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SUMMARY 

We  present  a  new  method  for  computing  the  impedance  matrix  (IM)  elements  in  the  method 
of  moments  for  geometries  described  by  bilinear  quadrilaterals  (BQ)  and  for  higher-order  basis 
functions.  Our  method  is  restricted  to  the  Electric  Field  Integral  Equation  and  focuses  on  the  self¬ 
elements  of  the  IM  and  elements  for  which  the  observation  point  is  near  the  integration  BQ.  The 
method  is  based  on  the  simple  idea  of  analytical  integration  along  one  of  the  BQ’s  parameters 
and  numerical  integration  along  the  remaining  one.  For  the  singular  (or  nearly  so)  part  of  the  in¬ 
tegral,  we  show  through  analysis  and  examples  that  our  method  can  provide  precision  to  1 5  sig¬ 
nificant  digits. 
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SECTION  1:  INTRODUCTION 

Electromagnetic  (EM)  modeling  and  simulation  (M&S)  is  increasingly  used  in  the  design 
stage  to  predict  the  behavior  of  electromagnetic  systems  in  diverse  and  complex  environments. 
Typical  examples  are  antenna  design,  in-situ  antenna-to-antenna  interaction,  radar  cross  section 
analysis  and  design,  EM  interference  and  vulnerability,  and  many  others. 

The  computational  electromagnetics  (CEM)  methods  that  are  used  in  developing  EM  M&S 
software  tools  divide  into  two  broad  categories:  those  that  operate  in  the  time  domain  and  those 
in  the  frequency  domain.  In  the  latter  case,  we  can  distinguish  again  two  categories:  methods  that 
approximate  the  physics  of  the  problem,  and  methods  that  do  not.  In  turn,  exact-physics  methods 
(as  the  latter  category  is  known)  are  split  into  two  categories:  those  that  operate  in  a  volume  set¬ 
ting  and  those  that  operate  over  surfaces.  The  latter  methods  are  known  as  boundary  integral 
equation  (BIE)  methods  and  the  present  study  falls  into  this  category.  The  numerical  solution  of 
the  BIE  is  usually  accomplished  in  CEM  using  the  method  of  moments  (MoM)  [1]  whereby  the 
integral  equation  is  converted  into  a  system  of  linear  equations  whose  unknowns  are  related  to 
the  currents  flowing  on  the  surface(s)  of  the  object  of  interest. 

In  MoM,  the  surface  of  the  object  of  interest  is  represented  by  a  number  of  geometric  pan¬ 
els.  Currents  are  defined  on  each  panel  and  the  total  current  on  the  object’s  surface  is  expressed 
as  a  linear  combination  of  the  panel  currents.  Thus,  in  practice,  the  actual  geometry  of  an  object 
(< e.g aircraft,  antenna,  terrain)  is  replaced  by  a  finite  number  of  surface  panels:  the  totality  of 
which  comes  close  to  the  true  geometry  of  the  object.  To-date,  most  MoM-based  EM  M&S 
commercial  tools  represent  the  actual  geometry  by  means  of  flat,  triangular  panels.  The  solution 
method  employed  after  the  triangular  discretization  is,  usually,  that  of  Rao,  Wilton,  and  Glisson 
[2].  The  drawback  of  using  flat,  triangular  panels  is  that  a  few  million  may  be  required  in  faith¬ 
fully  reproducing  a  large  and  complex  geometry.  Since  the  number  of  unknowns  in  the  resulting 
system  of  equations  is  related  to  the  number  of  triangle  edges,  this  may  result  in  a  system  with 
millions  of  unknowns,  requiring  substantial  computer  resources  and  long  execution  times.  It  is 
obvious  then  that  use  of  higher-order  panels  (panels  representable  by  higher-order  polynomials) 
could  result  in  a  large  reduction  in  the  number  of  panels  that  represent  the  actual  geometry.  This 
could  lead  to  a  substantial  reduction  in  the  size  of  the  system  of  equations.  A  very  well  written 
review  of  efforts  in  this  direction  up  to  2008  is  given  by  Notaros  [3].  At  this  time  (2015),  howev¬ 
er,  EM  M&S  software  based  on  the  methods  in  [3]  are  still  in  the  research  stage  with  one  excep¬ 
tion:  software  based  on  representing  geometry  using  bilinear  quadrilaterals  (BQ)  [4],  [5],  BQs 
are  surfaces  in  space  that  are  represented  in  terms  of  two  parameters.  The  representation  is  linear 
in  each  of  the  parameters  but  it  also  contains  their  product.  They  are  the  simplest  curved  panels 
in  use  in  CEM. 

An  important  part  of  the  application  of  MoM  to  an  object  comprising  a  number  of  BQs  is 
the  calculation  of  the  matrix  elements  of  the  resulting  system  of  linear  equations.  In  the  case  of 
the  electric  field  integral  equation  (EFIE)  ([1],  [2]),  the  matrix  element  consists  of  two  iterated 
integrals.  The  inner  integral  involves  the  product  of  the  free-space  Green’s  function  for  the 
Helmholtz  equation  multiplied  by  an  appropriate  basis  function  and  integrated  over  the  BQ  to 
which  the  source  point  belongs.  The  outer  integral  involves  the  result  of  the  first  integration, 
multiplied  by  a  test  function  and  integrated  over  the  BQ  to  which  the  test  function  belongs. 
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It  is  the  inner  integral  of  the  EFIE  matrix  element  that  we  study  in  this  report.  The  free- 
space  Green’s  function  is  of  the  form  cxp(ikR)  /  R,  where  k  is  the  wavenumber  while  R  is  the  dis¬ 
tance  between  the  source  point  on  the  integration  BQ  and  the  observation  point  (OP).  Thus,  the 
inner  integral  has  an  integrable  singularity  when  the  OP  lies  on  the  integration  BQ.  In  tenns  of 
numerical  integration  using  finite  precision,  this  singularity  manifests  itself  not  only  when  R  is 
equal  to  zero  but,  also,  when  it  is  close  to  it. 

There  is  a  considerable  amount  of  literature  on  the  evaluation  of  this  integral  [3],  [4],  [6].  In 
[6],  the  authors  consider  quadrilaterals  of  arbitrary  order  and  present  a  numerical-integration 
scheme  based  on  Duffy's  method  [7],  The  scheme  is  valid  only  if  the  OP  lies  on  the  integration 
QL.  They  then  proceed  to  compare  their  method  to  four  other  methods:  singularity  extraction 
([3],  [4]),  polar-coordinate  transfonnation  [8],  and  quadratic  [9]  and  cubic  [10]  rectangular  trans¬ 
formation. 

All  of  the  methods  just  mentioned  hold  for  arbitrary  degree  QL  panels.  In  the  present  study, 
we  focus  on  a  BQ.  This  allows  us  to  take  advantage  of  the  functional  form  of  the  distance  func¬ 
tion  and  develop  fonnulas  that  compute  the  singular  or  near-singular  part  to  machine  precision. 
In  Section  2,  we  review  a  BQ's  geometry.  In  Section  3,  we  present  three  QLs  that  we  use  to  test 
our  method.  In  Section  4,  we  present  arguments  why  some  of  the  present  approaches  may  not 
yield  the  precision  that  a  specific  problem  may  require.  In  Section  5,  we  examine  the  integral  un¬ 
der  consideration  and  we  split  it  into  two  integrals:  one  that  is  proper  (not  singular)  and  can  be 
evaluated  numerically  using  standard  methods,  and  an  improper  one  that  contains  an  integrable 
singularity.  We  then  show  that  the  latter  integral  (which  represents  a  Newtonian  potential  of  pol¬ 
ynomial  density  in  the  two  variables  of  integration)  can  be  reduced  through  repeated  integrations 
by  parts  with  respect  to  one  of  the  variables  to  a  double  integral  over  the  BQ  of  a  Newtonian  po¬ 
tential  of  density  one  and  a  series  of  single  integrals  that  are  not  singular  and,  hence,  can  be 
computed  through  standard  numerical  methods.  In  Section  6,  we  deal  with  the  remaining  singu¬ 
lar  integral  as  an  iterated  double  integral  in  the  BQ's  variables  and  we  show  that  we  can  evaluate 
either  of  the  two  single  integrals  analytically.  The  remaining  single  integral  contains  a  logarith¬ 
mic  singularity,  as  expected.  Through  numerical  examples  using  two  quite  distinct  quadratures, 
we  show  that,  as  the  OP  approaches  the  integration  BQ,  we  keep  losing  significant  digit(s)  (SD) 
for  the  remaining  integral.  At  this  point  it  appears  that  the  present  approach  does  not  work  as 
well  as  expected.  It  is  worth  mentioning,  however,  that  the  precision  we  obtain  does  not  go  be¬ 
low  eight  SD. 

In  Section  7,  we  revisit  the  result  of  Section  6  and,  upon  close  examination;  we  are  able  to 
pinpoint  the  reason  for  the  loss  of  SD.  The  logarithm's  argument  involves  the  difference  of  two 
numbers  that  are  almost  identical  at  the  origin.  This  leads  to  a  numerical  instability  and  eventual 
loss  of  SD.  By  manipulating  the  logarithmic  argument,  we  are  able  to  stabilize  it  and  show  that 
the  last  remaining  single  integral  can  be  evaluated  numerically  to  double  precision  (DP)  (15  SD) 
for  the  OP  not  only  close  but  on  the  BQ  itself.  In  Section  8,  we  use  the  third  BQ  with  additional 
OPs  on  it  to  provide  further  evidence  of  the  validity  of  our  approach.  In  Section  9,  we  deal  most¬ 
ly  with  irregular  BQs  (Section  2)  and  we  perform  sensitivity  studies  on  our  algorithms  to  deter¬ 
mine  how  close  we  can  get  to  an  irregular  BQ  before  we  start  losing  precision. 
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The  integration  routines  we  used  in  our  calculations  are  from  Mathematica®  7  [11]:  the 
Gauss-Kronrod  Quadrature  (GKQ)  and  the  Double-Exponential  Quadrature  (DEQ),  also  known 
as  the  Tanh-Sinh  Quadrature  (TSQ).  We  set  the  maximum  number  of  subdivisions  at  12  and 
asked  for  15  SD. 
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SECTION  2:  GEOMETRY  OF  A  BILINEAR  QUADRILATERAL 

In  this  section,  we  provide  a  basic  description  of  a  BQ.  With  respect  to  a  rectangular  coordi¬ 
nate  system  xyz,  we  define  a  BQ  by 

r (p,q)  =  ^[ru(l  -  P)( 1  ~q)  +  ri2(!  “  P)(.Q  +0  +  r2i (p  + 1)(1  ~q)  +  r22{p  + 1 ){q  +  D] , 

\p\<l,  \q\<  1.  (2.1) 

The  four  constant  vectors  r(/-  denote  the  four  comers  of  the  BQ.  We  note  that  (-1,  -1)  maps  to  rn, 
(-1,  1)  maps  to  ri2,  (1,  1)  maps  to  r22  and  (1,  -1)  maps  to  r2i.  We  show  this  in  Fig.  2.1. 


Figure  2.1.  An  example  of  how  the  basic  square  in  the  pq- plane  maps  to  a  BQ  in  the  xv-plane 
that  occupies  the  region  |x|  <  3 ,  |_y|  <  2 .  The  point  A  maps  to  1 ,  the  point  B  to  2,  point  C  to  3 
and  point  D  to  4. 

We  collect  terms  in  p,  q,  and  their  product  in  (1 . 1)  to  write 

r(  P,  q)  =  %)  +  rp  p  +  vqq  +  rpqPq ,  |p|<l,  \q\  <  1  (2.2) 

where 


r00=^[rll+r12+r21+r22] 

YP  =  ^[“rll  -  r12  +  r21  +  r22]  .  rq  =  ^[-rll  +  r12  “  r21  +  r22] 
rP?=^[rll-r12-r21+r22]- 


(2.3) 

(2.4) 

(2.5) 
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For  the  example  of  Figure  2.1:  r00  =  0 ,  r  =  3x ,  r  =  2 y  ,  r  =  0  . 

The  function  in  (2.2)  maps  points  from  the  square  \p\  <  1 ,  q  <  1  to  points  in  the  xyz- space. 
The  matrix  of  this  transfonnation  is 


M  = 


d  d  d 

-(.r-r)  — (y-r)  —  (z-r) 
op  op  op 

d  d  d 

-(xt)  — (y-r)  — (z-r) 
dq  dq  dq 


x-(rv+rpqq)  y-(rp+rpaq)  z-(rp+rpaq ) 


pep 


pep 


X-(ra+rpqP)  y-(ra+ruaP)  Z  ’  (ra  +  YvaP) 


(2.6) 


pqt 


pq1 


and  it  must  be  of  rank  2  if  we  are  not  to  have  any  singular  points.  This  means  that 

Mm)  x  8r(p,g)  Q 
dp  dq 


(2.7) 


at  every  non-singular  point  (p ,  q).  This  follows  from  the  fact  that  the  three  2x2  sub-matrices 
must  have  a  non-zero  detenninant.  Since  the  first  of  these  vectors  is  tangent  to  a  constant  g-curve 
while  the  second  to  a  constant  p- curve,  condition  (2.7)  says  that  the  two  vectors  should  not  be 
collinear  at  a  point  ([12],  pp.  56  -  57).  What  does  this  mean  geometrically?  From  (2.2)  we  have 
that 


dr(p,q)  dr(p,q) 

- X - 

dp  dq 


(rp  +  rpqCl)  X  (rq  +  rpqP)  =  rp  X  *q  +  V pq  X  (r^  "  ^ *pP)  • 


(2.8) 


If  we  assume  that  vp  and  vq  are  not  collinear  (and,  hence,  their  cross  product  is  not  zero),  then 
they  define  a  plane  in  space.  The  vector  (rqq  -  r pp^j  lies  on  that  plane.  Then,  a  necessary  condi¬ 
tion  for  (2.8)  to  be  equal  to  zero  is  that  the  vector  xpq  does  not  have  a  component  normal  to  that 
plane.  Differently  stated,  for  (2.8)  to  become  zero,  we  must  have  the  last  three  vectors  in  (2.2)  ly¬ 
ing  on  the  same  plane.  In  such  a  case,  (2.2)  describes  a  planar  BQ.  From  this  discussion,  it  is 

clear  that  the  test  for  determining  whether  a  BQ  is  planar  is  that  rpq  ■  (rp  x  rq  j  =  0  .  It  is  obvi¬ 
ous  from  this  condition  that  if  rp  and  rq  are  collinear,  we  then  have  an  irregular  BQ.  We  substan¬ 
tiate  these  statements  with  mathematics  in  Appendix  A. 

Although  (2.8)  can  be  equal  to  zero  only  for  a  planar  BQ,  the  converse  is  not  true,  i.e.,  (2.8) 
is  not  zero  for  all  planar  BQ.  We  investigate  this  further.  For  a  planar  BQ,  we  can  resolve  the  last 
vector  in  (2.2)  along  the  directions  of  the  preceding  two  vectors 

rpq  =  arp  +  Prq  (2.9) 


where  a  and  / 3  are  scalars.  In  place  of  (2.2),  we  now  have 

r(p,q)=roo+rpp(l  +  aq)  +  rqq(l  +  0p),  |p|<l,  |tf|<l. 


(2.10) 
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Using  this,  we  get 


dr(p,q)  ^  dr(p,q ) 

X 


dp  ^  =  [rp(\  +  aq)  +  Tqfiqy\rpap  +  rq(\  +  fip)] 

=  (r9  xrp)a/3Pq  +  (rp  xrq){l  +  aq)(\  + /3p)  =  (rp  xr^)(l  +  aq  +  /3p)- 


(2.11) 


If  the  cross  product  is  zero,  then  we  have  a  straight  line  in  space.  We  exclude  this  case.  When  is 
the  scalar  part  zero?  We  let 


P  =  Pp ,  Q-aq. 


(2.12) 


From  the  conditions  on p  and  q  in  (2.2),  we  have  that 

\P\  =  \P\\p\  —  \P\  >  |G|  =  |«||^|<|«|-  (2-13) 

We  also  set  the  scalar  factor  in  (2.1 1)  equal  to  zero 

l  +  Q  +  P  =  0 .  (2.14) 


We  want  to  know  whether  this  equation  is  satisfied  with  values  of  P  and  Q  in  the  range  of  (2.13). 
In  Figure  2.2,  the  lower  left  comer  of  the  rectangle  is  the  point  closest  to  the  straight  line.  The 
value  of  P  there  is  -  \P\  and  the  value  of  Q  for  the  straight  line  is  -  I  +  \p\.  For  no  intersection,  we 
want  -  1  +  \P  to  be  smaller  than  -  \a\;  for  one  intersection  point,  we  want  -  1  +  \P\  to  be  equal  to  - 
aj;  and,  for  an  uncountable  number  of  intersection  points,  we  want  -  I  +  \p\  to  be  greater  than  - 
\a\.  We  summarize  these  results  in  Table  2.1. 


Figure  2.2.  Geometry  for  the  intersection  of  the  straight  line  with  the  rectangle. 
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Table  2.1:  Solutions  of  (2.14)  subject  to  (2.13). 


t 

x  +  J3  <  1 

No  intersection  points 

No  solution  of  (2.14) 

Regular  BQ 

c 

x  +  f3  =  1 

One  intersection  point 

One  solution  of  (2.14) 

Irregular  BQ 

t 

x  +  J3  >  1 

Many  intersection  points 

Many  solutions  of  (2.14) 

Irregular  BQ 

We  demonstrate  these  results  using  three  BQs.  The  position  vector  for  the  first  is 

r(p,q)  =  rpp(l  +  0.2q)  +  rqq{l  +  0.3p),  (a  =  0.2 ,  (3  =  03)  (2.15) 

and  its  graph  is  shown  in  Figure  2.3.  This  is  a  regular  BQ.  The  position  vector  of  the  second  BQ 
is 

r(p,q)  =  rpp(l  +  0.5q)  +  rqq(l  +  0.5p)  ,  (a  =  0.5 ,  J3  =  0.5)  (2.16) 

and  its  graph  is  shown  in  Figure  2.4.  This  is  an  irregular  BQ  in  the  form  of  a  triangle.  The  posi¬ 
tion  vector  of  the  third  BQ  is 


r(p,q)  =  rpp(l  +  2q)  +  rqq{l  +  3p),  (a  =  2,j3  =  3).  (2.17) 

Its  graph  is  given  in  Figure  2.5.  This  is  also  an  irregular  BQ.  In  all  three  cases,  rp  =(3,0, -0.5) 
while  rq  =(0,2, -0.5)  . 


Figure  2.3:  Graph  of  (2.15);  a  planar,  regular  BQ. 
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Figure  2.4:  The  graph  of  (2.16);  a  planar,  barely  irregular  BQ. 


Figure  2.5:  The  graph  of  (2.17);  a  planar,  highly  irregular  BQ. 
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The  unit  normal  at  a  point  on  the  surface  is  given  by  ([12],  p.  62) 


8r(p,q)  dr  (p,q) 

X 


n  = 


dp 


dq 


dr(p,q)  dr(p,q)\ 

X 


(2.18) 


dp 


dq 


and  we  note  that  we  could  also  have  defined  it  as  the  negative  of  this.  The  surface  element  at  a 
point  is  defined  by 


clS  = 


dr(p,q )  dr(p,q ) 

X 


dp 


dq 


dpdq . 


(2.19) 


If  we  hold  one  of  the  parameters  constant,  we  can  measure  arc  length  along  the  other  ac¬ 
cording  to  the  formula  ([12],  p.  58) 


P  2 

=  1 


Pi 

Similarly, 


dr(p,q) 

dp 


Pi 


dP  =  j  \rp  +  rpqq\ dp  =  \rp  +  rpqq\ (p2  -Pl),  -\<pl<p2<\,  q  =  const. 

(2.20) 


Pi 


92 


9l 


dr(p,g) 

dq 


dq  -■ 


r  +  r 

q  pq 


p{<h-qi)' 


-l<ql<q2<l,  p  =  const. 


(2.21) 


We  point  out  that  if  we  have  a  curve  on  the  BQ  described  parametrically  by  p{t)  and  q(t),  then  we 
have  a  more  general  fonnula  for  its  arc  length  ([12],  p.  58).  In  any  case,  the  above  fonnulas  make 
sense  since  we  are  moving  along  straight  lines. 

We  consider  next  a  point  ( p0  q0  )  in  the  interior  of  the  square  of  definition  of  these  variables. 
This  maps  to  the  point  r0 ,  a  point  in  the  interior  of  the  BQ 

=  r  (  p() ,  q() )  =  r(X)  +  rp  p0  +  rqq0  +  rpq  p0q0 ,  |p0|<l,  |%|<1.  (2.22) 

We  can  reference  (2.2)  to  this  point  by  expanding  in  a  Taylor  series  about  it 

r(p,  q)=  r0  +  (rp  +  r pqq0  )(p~p0)  +  {rq  +  rpqPo  )(q-q0)  +  rpq(p- p0)(q-q0).  (2.23) 
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This  expression  is  exact.  We  can  easily  show  this  by  collecting  terms  in  powers  of  p,  q,  and  pq. 
At  the  point  r0,  we  define  coordinates  poqon.  The  unit  vectors  are 


dr(p0,q0) 

dr(p0,q0) 

dr(p0,q0) 

dr(p0,q()) 

fyo 

■ .  %  = 

dq0 

VI  —  - 

fyo 

X 

dq0 

dr(p0,q0) 

dr(p0,q0) 

?  11  — 

dr(p0,q0) 

„dr(p0,q0) 

fyo 

dq0 

fyo 

X 

dq0 

Although  Pq  and  %  are  perpendicular  to  n  ,  they  are  not  necessarily  perpendicular  to  each  oth¬ 
er.  Letting 


_ fr(P(h%) 

Po  “  a 

m 


„  _Sr(Po»«o) 

’  “  a 

dq0 


(2.25) 


we  can  write  for  (2.23) 

r ip,  q)  =  r0  +  p0w  +  Qov  +  rPquv  (2-26) 

where 


u  =  P~Pq,  v  =  q-q0;  -  (l  +  p0)  <  u  <  1  -  p0  ,  -  (l  +  q0)  <  v  <  1- q0.  (2.27) 

In  the  next  section,  we  will  describe  three  BQs  that  we  will  use  for  testing  purposes. 
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SECTION  3:  TEST  BILINEAR  QUADRILATERALS 

We  will  use  three  BQs  for  testing  our  ideas.  We  describe  them  below.  The  comers  of  the 
first  BQ  are  given  by 

in  =  (-3, -2,1)  ,  r22=(3,2,-l),  r12  =  (-3,2,0)  ,  r21=(3,-2,0)  (3.1) 

From  (2.3)  -  (2.5) 

r0  =  (0,0,0),  rp  =  ^(6,0, -l)  ,  r?  =|(0,4,-l)  ,  ^  =  (0,0,0).  (3.2) 

From  (2.2),  we  then  have 

r(p,q)  -  ~(6,0,-l)/7  +  ^(0,4,-l)g ,  |p|<l,  |r/  <  1 .  (3.3) 

The  graph  of  the  BQ  is  shown  in  Figure  3.1.  This  is  a  flat  BQ.  The  normal  is  given  by  (2.18). 


Figure  3.1:  The  BQ  of  Equation  (3.3). 

The  second  BQ  has  comer  vectors 

rn=  (-3-2,1),  ik  =  (3,2,1),  r12=(— 3,2,— 1),  r21  =  (3,-2,-l)  (3.4) 
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from  which 

r0=  (0,0,0),  rp  =  (3,0,0)  ,  rq  =  (0,2,0)  ,  r^=(0,0,l)  (3.5) 

so  that 

r(p,q)  =  (3,0,0) p  + (0,2,0) q  + (0,0,1) pq  ,  \p\  <  1 ,  <  1  -  (3.6) 

The  graph  of  this  BQ  is  shown  in  Figure  3.2. 


Figure  3.2:  The  BQ  of  Equation  (3.6). 

We  go  next  to  an  extreme  case 

rn=  (-3-2,10),  r22  =  (3,2,10),  r12  =  (-3, 2,-10) ,  r21  =  (3, -2,-10)  (3.7) 

so  that 

r0  =  (0,0,0),  rp  =  (3,0,0)  ,  rq  =  (0,2,0),  rp<?  =  (0,0,10)  (3.8) 

and 

r(/?, <7)  =  (3, 0,0)/? +  (0,2,0) <7 +  (0,0, 10)/*?,  |/?I<1,  |^|  <1.  (3.9) 

The  graph  of  this  BQ  is  shown  in  Figure  3.3. 
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We  note  that  condition  (2.7)  is  satisfied  at  all  points  of  these  QLs.  For  example,  for  the  last 
BQ  we  have 

x  =  K3,0,0) + (°’°’10W x  [(°'2’°)+ (°'0'10)  p] 

=  (3,0,0)  x  (0,2,0)  +  (3,0,0)  x  (0,0, 10)  p  +  (0,0, 10)  x  (0,2,0)^ 

=  6z-30py  -  20qx  (3.10) 

and  the  magnitude  of  this  vector  is  always  greater  than  zero.  This  means  that  we  have  a  tangent 
plane  at  every  point  of  the  BQ  with  a  normal  given  by  (2.18).  For  the  example 

,=  -20gx-30py  +  6z  (3  ,n) 

^(lOqf  +(30pf  +36 


-2  o  2 


Figure  3.3:  The  BQ  of  Equation  (3.9). 
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SECTION  4:  REVIEW  OF  EVALUATION  OF  AN  INTEGRAL 
We  justify  here  the  need  for  a  new  way  to  evaluate  the  integral 

1  1  its 

pim=\\p,(lm - dpdq  ,  l,m  =  0,1,2,...  (4.1) 

JiJ,  * 


over  a  BQ.  The  function  R  represents  the  distance  between  the  integration  point  (IP)  r  (on  the 
BQ)  and  another  point,  r',  that  we  call  the  OP.  This  integral  is  well  behaved  when  the  OP  is  some 
distance  away  from  the  BQ;  however,  when  the  OP  is  very  near  or  on  the  BQ,  then  the  integrand 
has  a  singularity  since  the  distance  function  R  can  become  equal  to  zero  or  very  nearly  so. 
Hence,  the  principal  issue  with  this  integral  is  how  to  calculate  it  when  the  OP  is  very  near  or  on 
the  BQ.  Two  principal  methods  exist  for  evaluating  (4.1).  One  relies  on  singularity  extraction 
and  the  use  of  an  auxiliary  tangent  plane  ([3],  [4])  while  the  other  on  coordinate  transfonnation 
([6],  [8],  [9],  [10]).  We  first  discuss  the  fonner  method. 

Kolundzija  and  Djordjevic  ([4],  p.  350,  2002)  advocate  the  following  splitting  of  the  integral 
so  as  to  isolate  its  singular  part 


p>»  =  J  J 


e!kR  | 

P'dm  —  ~  Po  V'  — 

K  K, 


1  1  1 

dpdq  +  p0'  qQm  J  J  —  dpdq,  l,  m  =  0, 1, 2, 


-i  -l  *o 


(4.2) 


where  ( p0,q0 )  corresponds  to  a  point  in  the  interior  of  the  BQ  and  is  the  projection  of  the  OP 

onto  the  BQ,  while  Ro  is  the  distance  between  the  OP  and  the  purely  linear  part  of  the  BQ’s  posi¬ 
tion  vector.  The  second  integral  can  be  evaluated  analytically.  For  the  first  integral,  the  authors 
write,  “The  integrand  of  the  first  integral... tends  to  zero  when  the  source  point  (IP)  approaches 
the  field  point  (OP).  Hence,  the  numerical  integration  of  this  integral  is  much  more  efficient  than 
the  numerical  evaluation  of  the  complete  ...  integral”.  We  will  show  below  that  this  statement 
does  not  stand  to  scrutiny. 

Notaros  ([3],  p.  2265)  mentions  the  same  splitting  of  the  integral  as  in  (4.2)  and,  for  the  first 
integral  in  (4.2)  suggests  “using  Gauss-Legendre  quadrature  formulas”.  What  is  most  likely 
meant  here  is  that  the  Gauss-Legendre  formulas  are  to  be  applied  along  the  p  and  q  coordinates. 

We  proceed  to  examine  these  statements  closely.  The  worst  case  scenario  is  when  the  OP  is 
on  the  BQ,  i.e.,  the  point  ( p0,q0 )  •  For  the  position  vector  to  any  point  on  the  BQ,  we  have  from 
(2.23)  that 


r ip,  q)  =  r  *  +  (r;,  +  r pqq0  )  {p  ~  P0  )  +  (r,  +  rpqPo )  {q  ~  q0  )  +  rpq  (p  -  p0 )  (q  -  q0  ) ,  r*  =  r  (p0,  q0  )  . 

(4.3) 
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We  let 


dr(p0,q0 )  ,  ^r(p0,g0) 

P(A>>tfo)  = - - - =  P0>  q(  A) ’</(.)  =  — 7 —  =  q0 


(4.4) 


<7Po 

and  we  also  make  the  change  of  variables  (2.27) 

u  =  p-  p0  ,  v  =  q-q0;  -(l  + p0)<u<l- p0  ,  ~(l  +  q0)<v  <l-q0.  (4.5) 

We  can  then  write  for  (4.3) 

r (P,  q)  =  r*  +p 0u  +  q0v  +  rpquv  .  (4.6) 

The  linear  part  of  this 

r0(p,q)  =  r*+p0u  +  q0v  (4.7) 

is  a  vector  to  a  point  on  the  plane  tangent  to  the  BQ  at  r*.  With  this  notation,  we  have  that 

/?  =  |r-r*| ,  R0  =|r0  —  r*| .  (4.8) 


For  the  first  integrand  in  (4.2),  we  write 
f{u,v)  =  (u  +  p0)'(v  +  q0) 


ikR(u.v)  i 

Po%m 


R(u,v)  1 0J0  R0(u,v ) 
and  proceed  to  evaluate  the  limit  at  the  origin.  We  let  first  u  go  to  zero 

1 


(4.9) 


/  (0,v)  =  p0‘  ( v  +  q0)" 


R(  0,v) 

From  (4.7)  -  (4.8),  we  get  that 

PikM 

f(0,v)  =  p0‘(v+qoy 


JkR(  0,v) 

c  l  m 

PoVo 


R0  (o> v) 


(4.10) 


Kv 


■  Po'%m 


Kv 


(4.11) 


We  proceed  to  expand  the  binomial 


/&M|q0|  m  . 

,z 

( in' 

V  Qo  |  n= 0 

w 

n  „  m—n  l  _  m 

v  7o  ~p0q0 


l 


pQ'q o'”  ■ 


AMI«iol 


-1 


AMkl 


+  Po 


-Z 


We  take  now  the  limit  as  v  goes  to  zero 


'IK 

vnq0 


n  m—n 


(4.12) 
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/  m— 1 

]imf(0,v)  =  ikp0'q0m±mp°q°  .  (4.13) 

v— >o  n 

|4o| 

where  the  plus  sign  corresponds  to  positive  values  of  v  and  the  minus  to  negative.  We  see  that 
the  imaginary  part  of  the  integrand  does  possess  a  limit.  This  is  to  be  expected  because  the  imag¬ 
inary  part  has  no  singularity  since  it  is  of  the  form  sin(x)/x.  The  real  part,  however,  does  not  pos¬ 
sess  a  limit  since  the  limit  value  depends  on  the  direction  from  which  the  point  in  question  is  ap¬ 
proached. 


We  next  let  v  be  equal  to  zero  in  (4.9) 

JkR(u,  0) 


f(u,0)  =  (u  +  p0)  qQ" 

=  (u  +  p0)'  q0m 


■  Po< 


R(u,0)  °  °  R0(u, 0) 


A|p(po.9o)|H 


P(Po^0)| 

e*|Po||“l  J  f  l  ^ 


„  l  „  m 

'  Po  clo 


1 


I 


Po  "  «jV“ 


n  l—n  l  _  m 

u  Po  ~  Po  9o 


|Po||“| 

i 


Po  M 


/  m 

_  Po  Ho  r^inbiH 


Po  M 


.  j~e'*|Poll«l  _  J  J  . 


I«*|p(Pb.9b)|M  / 


% 


z 


Po  «  »=1 


u"Po  H 


We  now  take  the  limit  as  u  tends  to  zero 


(4.14) 


i  l— l  m 

lim /  (u, 0)  =  ikp0'qQ'n  ±  Po.  q.°  (4.15) 

I  Pol 

with  the  dual  sign  as  in  the  previous  case.  Again,  the  imaginary  part  of  the  integrand  has  a  limit 
while  the  real  part  does  not. 


We  see  that,  for  the  real  part  of  the  integrand,  the  two  operations  produce  different  limits; 
moreover,  within  each  operation  the  limit  is  direction  dependent.  This  means  that  the  limit  of  the 
real  part  of  the  integrand  at  ( p0,q0 )  does  not  exist  and,  hence,  the  function  is  not  and  cannot  be¬ 
come  continuous  there.  In  terms  of  a  numerical  evaluation  of  the  first  integral  in  (4.2),  this  means 
that  we  need  a  superior  cubature  in  order  to  produce  good  results.  We  also  observe  that,  if  /  and 
m  are  equal  to  zero,  then  the  limit  does  exist  for  the  real  part  and  is  equal  to  zero.  Thus,  the  claim 
of  Kolundzija  and  Djordjevic  [4]  that  the  integrand  goes  to  zero  at  this  point  is  correct  only  for 
this  case.  We  present  additional  evidence  below,  based  on  numerical  calculations. 


In  performing  numerical  demonstrations,  we  use  the  three  BQs  we  presented  in  Section  3.  In 
all  three  cases,  we  place  the  singular  point  at 

po=qo  =  \/4  (4.16) 


so  that 
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u  =  p  — —  ,  v  =  q-—;  -1.25  <u<  0.75  ,  - 1.25  <  v  <  0.75  .  (4.17) 

4  4 

We  will  deal  only  with  the  real  part  of  the  integrand  since  the  imaginary  does  not  present  a  prob¬ 
lem. 


Re{fj(u,v)}  = 


c  n 

/ 

r  n 

l  4  J 

l  4 ) 

cos  kRj  (m,v)J 
Rj(u,v) 


'n 

Ay 


l+m 


R0J  (+V) 


j  =  1,2,3 


(4.18) 


with  the  subscript  j  referring  to  one  of  the  three  BQs.  Moreover,  from  the  expansion  of  the  cosine 
about  zero,  we  see  that  only  the  leading  term  is  singular,  while  the  rest  go  to  zero  with  Rj.  We 
dispense  with  those  tenns  and  write 


hj(u,v)  = 


(  0 

l 

(  0 

m  ^ 

rn 

u  +  — 

V  H - 

— 

— 

V  4  J 

v  4  J 

R:(u,V ) 

U  J 

l+m 


Roj(u’v)' 


j  = 1,2, 3. 


(4.19) 


This  is  the  surface  we  will  plot  about  the  origin. 

We  begin  with  the  flat  BQ  of  Figure  3.1.  The  position  vector  is  given  by  (3.3)  and,  hence, 


Ri(u,v) 


i  (6, 0, -1)  ii +  i(0, 4, -l)i 


:^|i6w  +  y4v-£(w-i-v)|  -  +  (4v)"  +  (n  +  v)" 

(4.20) 


Moreover, 


RiU  ( u,v)  =  Rl  (u,v). 


(4.21) 


From  (4.19) 


t\  (u,v) 


(  O' 

l  4  J 


(  0 

V  H - 

V  4  J 


f 

V 


n 

4y 


l+m 


_ \_ 

Rl  (u,v) 


(4.22) 


This  expression  is  equal  to  zero  when  /  =  m  =  0,  as  it  should.  In  Figs.  4. 1  and  4.2,  we  display  the 
cases  /  =  1,777  =  0  and  /  =  0,777  =  1  ,  respectively.  In  the  first  case,  we  observe  a  signum-like  be¬ 
havior  near  the  origin  and  along  the  w-axis  while,  in  the  second  case,  this  behavior  is  along  the  v- 
axis.  In  Figure  4.3,  we  display  the  case  /  =  1,77?  =  1  where,  as  expected  the  signum-like  behavior 
manifests  itself  along  the  bisector  of  the  two  axes  as  well  as  along  the  axes  themselves.  All  three 
figures  support  the  analysis  presented  above  and  the  surfaces  depicted  there  are  not  easy  to  du¬ 
plicate  numerically,  even  in  this  simple  case  of  a  flat  BQ.  Admittedly,  the  singularity  can  be  re¬ 
moved  through  a  rotation  of  the  BQ  so  that  the  normal  to  the  BQ  points  along  the  z-axis.  Once 
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this  is  done,  we  can  integrate  with  respect  to  one  of  the  variables  to  remove  the  singularity.  This 
can  be  done  only  for  a  flat  BQ. 


Figure  4.1:  Graph  of  the  surface  (4.22)  for  I  =  l,m  =  0 .  The  surface  has  a  signum-like  behavior 
near  the  origin  and  along  the  u-axis. 


Figure  4.2:  Graph  of  the  surface  (4.22)  for  I  =0,m  =  l.  The  surface  has  a  signum-like  behavior 
near  the  origin  and  along  the  v-axis. 
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Figure  4.3:  Graph  of  the  surface  (4.22)  for  /  =  l,m  =  1.  The  signum-like  behavior  near  the  origin 
manifests  itself  along  the  two  axes  and  along  a  45 -deg  line  with  respect  to  the  two  axes,  as  it 
should. 


We  proceed  to  check  the  values  in  these  graphs.  For  the  limits  in  (4. 1 1)  and  (4.13),  we  write 


lim/z,  (0,  v) 


m 

m 

l+m— 1 

lim h  (u, 0)  =  ±t — 7 

■  1  ^ 

|q0| 

v4y 

S 

0 

0 

(4.23) 


Using  (4.4)  and  (3.3),  we  write  for  these  limits 
2m 


lim/ij  (0,v)  =  ±-= 
v-*0  V  '  *Jl7 


1 

v4y 


=  ±0.485  lm 


1 

v4y 


(4.24) 


and 


lim/?.  (u,0)  =  ± 

u — >0  1  v  ' 


21 


f  1  V+m_1  f  i  A 


l+m-l 


Jyj 


=  ±0.3288/ 


(4.25) 


v^y 


For  the  case  /  =  1,777  =  0,  we  see  from  (4.22)  that  (4.24)  is  satisfied.  Moreover,  from  (4.22)  and 
(4.20) 


\  (u, 0)  = 


\(  0 

rivi 

— 

— 

LI  4 J 

U  )_ 

1 


2  u 


R{(u, 0)  y/37  \u\ 


(4.26) 
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which  is  in  agreement  with  (4.25).We  display  the  graph  of  this  function  in  Figure  4.4. 
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Figure  4.4:  Graph  of  the  function  /?,  (n,0) ,  as  defined  in  (4.22),  for  /  =  l,m  =  0 . 


Similar  comments  can  be  made  about  the  cases  /  =  0,m  =  1.  In  the  case  /  =  l,m  =  1,  we  get 
results  as  for  the  other  two  cases  but  we  have  one  more,  interesting  result  along  the  bisector.  Let¬ 
ting  v  =  u  in  (4.22),  we  get 


/ij  (ii,u)  = 


yfU 


,  1  U 

u  +  —  — 
2  \u\ 


u^O. 


(4.27) 


We  show  the  graph  of  this  in  Figure  4.5.  This  discontinuity  is  at  least  as  nasty  as  the  previous 
ones.  We  encounter  the  same  behavior  along  u  for  /  =  2,m  -  0  and  along  v  for  /  =  0,m  =  2 . 


Figure  4.5:  Graph  of  /?,  ( u,u )  for  /  =  l,m  -  1 . 
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So  as  not  to  crowd  this  section,  we  analyze  the  other  two  BQs  of  Section  3  in  Appendix  B. 
From  the  results  there  and  the  present  analysis,  we  conclude  that,  when  the  OP  is  near  or  on  the 
integration  BQ,  the  singularity  extraction  technique  is  not  sufficient  if  the  first  integral  in  (4.2)  is 
to  be  computed  in  pq-  or  wv-coordinates.  The  one  and  only  exception  is  when  both  /  and  in  are 
equal  to  zero. 

We  next  discuss  the  second  method,  that  of  a  coordinate  transfonnation.  For  simplicity,  we 
concentrate  on  the  transfonnation  to  polar  coordinates  ([3],  [8]).  If  the  OP  lies  on  the  integration 
BQ,  this  approach  works  very  well.  When,  however,  the  OP  is  removed  from  the  BQ,  it  becomes 
ineffective.  In  order  to  demonstrate  our  point,  we  use  the  second  BQ  of  Section  3  and  the  point 
of  (4.16).  Using  this  point,  we  erect  the  normal  along  which  we  designate  a  distance  h  from  the 
BQ.  This  is  where  the  OP  is  located.  We  also  introduce  polar  coordinates 

u-pcoscp,  v  =  psiri(p;  p>  0,  0<cp<27i  (4.28) 

with  u  and  v  defined  in  (4.5).  In  Figs.  4.6.1  -  4.6.4,  we  plot  the  ratio  pi  R  where  R  is  the  distance 
between  the  OP  and  the  source  (integration)  point  on  the  BQ.  From  these  graphs,  we  see  that  the 
integrand  is  not  very  smooth  near  the  origin.  By  switching  to  polar  coordinates,  we  remove  the 
(near-)  singularity  at  the  origin  but  we  are  still  left  with  a  function  that  is  not  easy  to  integrate  to 
high  accuracy. 

To  better  grasp  the  behavior  near  the  origin,  we  have  also  taken  cuts  along  u  =  0  and  v  =  0. 
We  exhibit  these  graphs  in  Figs.  4.7.1  -  4.7.5.  From  the  first  four  graphs,  we  see  that,  for  h  dif¬ 
ferent  from  zero,  we  have  continuity  at  the  origin  and  the  value  there  is  zero.  When  h  =  0,  how¬ 
ever,  the  limit  at  the  origin  does  not  exist,  as  we  see  from  Figure  4.7.5.  We  see  also  that,  the  larg¬ 
er  the  h  is  (in  absolute  value),  the  more  we  have  to  worry  about  the  accuracy  of  the  result.  As  h 
gets  small,  the  gap  between  the  two  sides  of  the  curve  diminishes  dramatically,  making  the  nu¬ 
merical  evaluation  easier.  All  in  all,  we  conclude  from  this  brief  discussion  that,  even  when  we 
use  polar  coordinates,  we  have  to  be  careful  how  we  develop  numerical  integration  algorithms. 
We  note  that  an  analytical  evaluation  cannot  be  perfonned  on  either  the  radial  or  the  angular  in¬ 
tegral. 
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Figure  4.6.1.  Graph  of  the  integrand  pi  R  for  h  =  0.001  about  the  origin  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16).  To  better  illustrate  its 
properties,  we  show  the  graph  in  an  upside-down  position. 


Figure  4.6.2.  Graph  of  the  integrand  pi  R  for  h  =  10'6  about  the  origin  for  the  sec¬ 
ond  BQ  of  Section  3.  The  origin  is  the  point  (4.16).  If  we  turn  this  graph  upside 
down,  we  get  a  figure  that  resembles  Fig.  4.6.1.  Conversely,  if  we  turn  Fig.  4.6.1 
upside  down,  we  get  a  figure  that  resembles  the  present  one. 
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12 

Figure  4.6.3.  Graph  of  the  integrand  p  /  R  for  h  —  10'  about  the  origin  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 


Figure  4.6.4.  Graph  of  the  integrand  p  /  R  for  h  =  0.0  about  the  origin  for  the  sec¬ 
ond  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 
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Figure  4.7.1.  Graph  of  the  integrand  p  /  R  for  h  =  0.1  along  v  =  0  (left)  and  it  =  0  (right)  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 


Figure  4.7.2.  Graph  of  the  integrand  p  /  R  for  h  =  0.1  along  v  =  0  (left)  and  u  =  0  (right)  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 


Figure  4.7.3.  Graph  of  the  integrand  pi  R  for  h .  =  10'6  along  v  =  0  (left)  and  u  =  0  (right)  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 
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12 

Figure  4.7.4.  Graph  of  the  integrand  pi  R  for  h  =  10'  along  v  =  0  (left)  and  u  =  0  (right)  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 
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Figure  4.7.5.  Graph  of  the  integrand  p  /  R  for  h  =  0.0  along  v  =  0  (left)  and  u  =  0  (right)  for  the 
second  BQ  of  Section  3.  The  origin  is  the  point  (4.16). 


The  next  question  we  may  ask  is:  what  if  we  use  polar  coordinates  and  the  tangent -plane  ap¬ 
proximation?  We  employ  (4.8)  and  compute  the  integrand  p  ((i  /  r)  -  (l  /  r0  ))  .  We  display  the 

graph  of  this  function  in  Figs.  4.8.1  -  4.8.5.  We  see  that  the  graphs  have  no  singularities  and  that 
they  are  relatively  smooth;  in  fact,  it  appears  that  the  smaller  the  h,  the  smoother  the  graph.  Since 
the  integral  of  1  /  Ro  can  be  evaluated  analytically,  this  approach  may  lead  to  good  accuracies  in 
computing  (4.1).  In  the  following  sections,  we  propose  an  approach  in  rectangular  coordinates 
that  makes  the  use  of  the  auxiliary  tangent  plane  unnecessary. 
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Figure  4.8.1.  The  integrand  p((l//?)-(l//?0))  for  the  second  BQ  of  Section  3 
and  for  h  =  0. 1 .  The  origin  is  the  point  (4.16). 


Figure  4.8.2.  The  integrand  p((l//?)-(l//?0))  for  the  second  BQ  of  Section  3 
and  for  h  =  0.001.  The  origin  is  the  point  (4.16). 


29 


NAWCADPAX/TR-20 15/214 


Figure  4.8.3.  The  integrand  /?((l//?)-(l//?0))  for  the  second  BQ  of  Section  3 
and  for  h  =  10'6.  The  origin  is  the  point  (4.16). 


Figure  4.8.4.  The  integrand  p((l//?)-(l//?0))  for  the  second  BQ  of  Section  3 

12 

and  for  h  =  10'  .  The  origin  is  the  point  (4.16). 
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Figure  4.8.5.  The  integrand  p((l//?)-(l//?0))  for  the  second  BQ  of  Section  3 
and  for  h  =  0.0.  The  origin  is  the  point  (4.16). 
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SECTION  5:  EXAMINATION  OF  ORIGINAL  INTEGRAL 

In  this  section  we  take  a  closer  look  at  the  integral  in  (4.1).  We  show  that  we  can  split  the  in¬ 
tegral  into  two  parts:  one  with  an  analytical  integrand  and  one  with  an  integrable  singularity.  For 
the  latter,  we  show  that  we  can  employ  integration  by  parts  to  perform  analytically  the  integra¬ 
tion  with  respect  to  one  of  its  variables,  leaving  an  integral  with  respect  to  the  other  variable 
whose  integrand  contains  at  worst  a  logarithmic  singularity.  In  subsequent  sections,  we  will 
show  that  the  remaining  integral  is  computable  to  machine  precision.  We  begin  by  re-writing 
(4.1)  in  the  form 

Plm  —  J  J  p  qm - ^  ^  dpdq  +  ik  J  J  p'qm  ^  - dpdq  ,  l,m  =  0,1,2,...  (5.1) 


Since  the  function  sin(.v)/.r  can  be  computed  to  all  available  precision  for  any  value  of  x,  no  mat¬ 
ter  how  small,  the  integral  of  the  imaginary  part  can  be  computed  using  established  quadratures. 
Most  of  the  contribution  will  come  from  the  neighborhood  of  R  =  0.  We  dispense  then  with  the 
imaginary  part  of  (5.1)  and  concentrate  on  the  real  part 


R e{iU=JJpY! 


cos 


(«) 


R 


dpdq,  l,  m  =  0,1,2,... 


(5.2) 


The  integrand  has  an  integrable  singularity  when  R  is  equal  to  zero  or,  for  finite-precision  arith¬ 
metic,  in  a  neighborhood  of  it.  Just  as  in  Section  4,  we  assume  that  the  OP  is  an  interior  point  (po , 
<7o)  of  the  BQ  and  employ  the  notation  of  (4.3)  through  (4.8).  As  we  shall  show  below,  the  point 
can  be  anywhere  on  the  BQ,  including  the  BQ’s  boundary. 

We  proceed  to  work  on  (5.2)  and  write 


Re(/M  = 


=  J  J  Pd"  COS^  ^  dpdq  +  J  J 


11  In 

p  q 


R 


dpdq 


=  -2  J  J  p'qm  sin  (kR/2)  dpdq  +  J  J  P_^_dpdq 


-1-1 
1  1 


R  JiJi  R 


J  J  pqm  (kR) 


-1-1 


sin  (kR/  2) 
kR/  2 


dpdq  +|| 


11  In 

p  q 


-\  -i 


R 


dpdq 


(5.3) 


and  substitute  the  result  in  (5.1) 


i  i 


Pim=--\\pqm  (kR) 


-1-1 


sin(k/?/2) 


=  ik\  J  P  q 


kR/ 2 

sin  (kR  /  2)  cos  (kR  /  2) 


1  ’A"  Lo«sinW 


dpdq  +||  —  dpdq  +  ikj  J  plq 


-i-i 


-i-i 


kR/ 2 


k  1  1 

dpdq 11  pl  q"  (kR) 


-1  -1 


-1-1 
“|2 


sin(kR/2) 
kR/ 2 


kR 
1 

dpdq  +  J  J 


dpdq 


11  In 

p  q 


-i-i 


R 


dpdq 
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=  ik\  j  p'q 


sin  (kR  /  2) 
kR/2 


11  /  .. 

p  q 


[cos  (kR  /  2)  +  i  sin  (kR/2  )  J  dpdq  +  II  — — — dpdq  (5.4) 


or 


Pbn  =ik\  J  pq 


sin  (kR  /  2) 
kR/2 


ii  / 

p  q 


e,kRldpdq  +  JJ — — — dpdq,  l,m  —  0,1,2,...-  (5.5) 


The  first  of  these  integrals  has  no  singularity  and  can  be  evaluated  numerically  using  stand¬ 
ard  methods.  We  note  that  the  since  function  emphasizes  the  values  of  the  integrand  about  R  =  0 
and  de-emphasizes  those  away  from  it;  moreover,  the  frequency  has  been  cut  in  half,  resulting  in 
a  less  oscillatory  integrand. 

To  work  on  the  last  integral,  we  need  the  expression  for  the  distance  function.  We  will  as¬ 
sume  that  the  OP  projects  to  a  point  (p0,q0)  of  the  BQ.  From  (2.23),  the  position  vector  to  a 

point  ( p0 ,q0  )  of  BQ  is 

r (P,  q)=  r0  +  (r„  +  rpqq0 ) ( P  ~  PQ ) +  (r,  + : rpqPo )  (q- -  q0 ) + : rpq  (p  -  p0 ) (q  -  q0 ) .  (5.6) 

Using  (2.25),  we  can  write  this  in  the  form 

r  O.  q)  =  ro  +  Po  (p  -Po)  +  q<)  -*2o)  +  vPq  (p-Po)(4-4o)-  (5-7) 

This  is  the  IP  on  the  BQ.  The  vectors  po  and  qo  define  a  plane  tangent  to  the  BQ  at  the  origin.  We 
define  the  OP  by 

r'  =  r0  +  Po  Xq°  h  =  r0  +  nh  .  (5.8) 

Po x  q(, 


With  (2.27),  we  can  then  write  for  the  distance  function, 

R (u,  v)~  =  |r  - r'[  =  r2  -  2r  •  r'  +  r 2  =  |p0w  +  q0v  +  r/J(/wv|  -  2 (p0M  +  q0v  +  ^pquv)  ■  nh  +  h 1 

=  |Pow  +  <lovr  +  M  (MV)2  +  2(Pow  +  qov)  •  rPquv  -  2(Pow  +  %v  +  rPquv)  ■  nh  +  h2  .(5.9) 
We  return  to  the  last  integral  in  (5.5)  and  write  for  it 


1  1  l  m 

p  q 


I  (p0,q0,l,m)=  f  f - dpdq,  l,m  =  0,1,2,... 

j  j  p 


(5.10) 


With  the  transformation  (2.27),  we  get  integrals  of  the  form 


l-<70  1  ~P0 

l(p0,q0,l,m)=  J  J  ——dudv , 

-(1+9o)-(1+ro) 


l,m  =  0,1,2,... 


(5.11) 
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We  show  below  how  this  double  integral  can  be  reduced  to  a  single  integral. 
Case  1:  /  =  m  =  0. 


r/  \  ‘7°  '7°  dudv 

7(p0,<?0,0,0)  =  7  =  J  J  — .  (5.12) 

“(l+9o)  “(l+p0) 

The  distance  function,  given  by  (5.9),  can  be  considered  as  a  second  degree  polynomial  in  u  or  in 
v.  As  such,  one  of  the  two  integrals  can  be  evaluated  analytically  ([13],  p.  81,  Formula  2.261). 
The  resulting  expression  involves  the  logarithm  of  a  complicated  function  of  the  remaining  vari¬ 
able.  This  function  goes  to  infinity  when  the  remaining  variable  is  zero  and  as  h  tends  to  zero. 
We  shall  study  this  integral  in  detail  in  subsequent  sections.  As  an  example,  we  consider  the  dis¬ 
tance  as  a  function  of  u.  From  (5.9),  we  write 

/?(m,v)  =  ^ jc 2 ( v)z/ 2  +  7?(v,/z)m  +  A(v,/z)  ,  C(v)>0  (5.13) 


where 


A(v,h)  =  \q0\2  v2  +h2 ,  B(v,h)  =  2[q0 -(p0  +rpgv)-n-rpqh]v  ,  C(v)  =  |p0  +rp?v| .  (5.14) 


If  we  integrate  (5.12)  with  respect  to  u,  we  get 


!-?o 

/=  j  In 

“(1+9o) 


2C  (v)  7?  (1  -  p0 ,  v)  +  2C2  (v)  (1  -  pQ )  +  B  (v,  h  ) 

dv 

2C(v)7?(-(l  +  /r0),v)-2C2(v)(l  +  /?0)  +  7?(v,/7) 

C(v) 

(5.15) 


We  can  easily  verify  that  the  denominator  of  the  logarithmic  argument  becomes  zero  at  v  =  0  and 
h  =  0.  The  numerator  becomes  zero  when  pQ  =  1  and  v  and  h  are  both  zero. 


Case  2:  /  =  0,  m  =  1. 


1  la  1  Pa 

7(p0,r/0,0,1)=  |  |  —  (5.16) 

~(1+9o)  -(1+Po) 

In  this  case,  we  proceed  as  in  Case  1  to  get  the  same  expression  as  in  (5.15)  but  with  the  inte¬ 
grand  multiplied  by  v.  The  limit  for  the  integrand  as  v  tends  to  zero  and  h  =  0  now  exists  and  is 
zero.  Thus,  in  the  numerical  evaluation  of  the  outer  integral  in  (5.16),  all  we  need  do  is  to  define 
the  value  of  the  integrand  at  zero  to  be  zero.  If  /  =  0  and  m  =  1,  then  we  consider  the  integral  with 
respect  to  v  as  the  inner  integral  and  proceed  as  above. 
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Case  3:  l  >  1,  m>  / . 
The  integral  is 


For 


l-?o  l~Po  l 

I(p0,q0,l,m)=  [  dvv”'  [  =du 

-(4)  -(4)  ijc  (v)u  +B (v, h ) M  +  A (v, h ) 

1  '"j?0  vm  1_f°  m,_1  [" 2C2  (v)zr  +  5 

=  -  dv  ,  L  J  -.du 

2  -(4)  C  (v)  -(4o)  a/C2  (v)«2  +  fi(v,/z)w  +  A(v,h) 

1  '7°  vmB(v,h)  ‘7°  7  1 

—  dv - A  ■ 

2-(i+9o)  C  (v)  -(Up^yjc2  (v)u2  +B(v,h)u  +  A(v,h) 

the  first  inner  integral  on  the  right-hand  side,  we  can  write 

l-f  ,/-'[2C=(v)a  +  B(v,A)] 

MA»tfo>0=~  J  ,  =  <*“=“  S(H-V) 

7  )^C  (v)w  +B(v,h)u  + A(v,h) 


du  . 


(  l+Po)  ' 


»=l-Po 

M=-(1+/)0) 


i-Po  , _ 

-(/-l)  |  ul~2 ^Jc2  (v)u2  +  B(v,h)u  +  A(v,h)du  . 

“(l+Po) 


(5.17) 


(5.18) 


The  last  integral  in  this  expression  can  be  evaluated  using  the  recursion  formulas  in  Gradshteyn 
and  Ryzhik  ([13],  Section  2.26,  p.  81).  From  ([13],  Formula  2.260.1,  p.  81),  we  have  that 


l-Po 


u'~3R 3 


[  u'  2R(u,v)du  ,  ,  , 

-L  lcld) 


(2f-3)B(v,ft) 

»=-(l+p0) 


1-Po 


9  .  .  f  u1  3R(u,v)du 

2lc-(d  -L 


T  u“R(u,v)du. 

<L) 


(5.19) 


If  we  use  this  formula  repeatedly,  we  end  up  with  a  number  of  rational  functions  of  v,  with  no 
poles  in  the  interval  of  integration,  and  an  integral  of  the  distance  function  multiplied  by  a  ra¬ 
tional  function  of  v.  From  ([13],  Formula  2.262.1,  p.  82),  we  have  that 


— 


R(u’v)  4(4) 


4  A{v,h)C2{v)-B2{v,h) 
8 C2  (v) 


1— Po 


du 


(5.20) 


We  will  show  how  to  compute  the  last  integral  in  the  next  section. 
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For  the  inner  integral  of  the  second  integral  on  the  right-hand  side  of  (5.16),  we  see  that  we 
have  an  integrand  of  exactly  the  same  form  as  the  integral  we  started  with,  except  that  the  power 
of  u  has  been  reduced  by  one.  We  can  then  apply  the  same  rule  to  this  integral. 

If  m  >  1  and  /  >  m,  then  we  consider  the  integral  with  respect  to  v  as  the  inner  integral  and 
proceed  as  above.  In  conclusion,  we  see  that  the  integral  we  must  study  is  the  integral  I(p 0,  qo, 
0,0)  in  (5.12).  We  proceed  to  do  that  in  the  next  section. 
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SECTION  6:  EVALUATION  OF  THE  INTEGRAL  OF  HR  IN 
RECTANGULAR  COORDINATES 

As  we  demonstrated  in  the  last  section,  the  integral  in  (5.11)  can  be  reduced  to  a  number  of 
well-behaved  single  integrals  and  an  integral  of  the  form 


!-«0  i-Po 


1 dv  1 


du 


-l-9b 


(6.1) 


We  have  omitted  from  the  integrand  a  function  of  v  that  is  at  least  continuous  since  its  presence 
or  absence  does  not  affect  our  argument.  We  examine  here  how  well  we  can  compute  this  inte¬ 
gral.  From  (5.7)  and  (5.8),  we  have  that 


R («,  v)  =  |r  Ar'| ,  r  =  r0  +  p 0u  +  q0v  +  r  uv ,  r'  =  r0  + 


,  PoXq() 


Poxq(» 


h  =  r0  +  nh  . 


(6.2) 


The  definitions  of  the  various  symbols  are  as  in  (2.23)  and  the  equations  following  it.  The  IP  is  r 
while  the  OP  is  r'.  The  OP  is  near  the  BQ  and  projects  to  a  point  of  the  BQ  that  corresponds  to 
the  point  (po,  go).  This  point  is  the  new  origin  of  rectangular  coordinates  (u,  v). 


We  consider  now  the  inner  integral 


l-Po 


du 


(6.3) 


For  the  distance  function,  we  write 


R2  (u,v)  =  |r|_  -  2r  •  r'  +  |r'|"  =  |p0zr  +  q0v|"  +  2(p0w  +  q0v)  •  r pquv  +  (uv)  - 2rpq  ■  nhuv  +  h 

=  |po|2  «2 + 2Po  •  q0wv + |q0|2  v2 + 2(po» + q0v  -”h)  •  rPquv + \rPq\2  (uvf + h 2 

=  (|Po|2  +  2Pn  ■  V  +  \rpq\  V2) “2  +  2[Po  •  q0  +  (q0v  -  nh)  ■  rpq ] vu  +  |q0|~  v2  +  h2 

=  |p0  +  W’\2  u 2  +  2 [q0  •  (po  +  rPqv )  -  n  ■  rPqh\ vu  +  |q0|2  V2  +  h2  (6-4) 


or 


R(u’v)  =  J|Po  +  vpqv[ 1,2  +  2  q<)  ■  (Po  +  rPqv ) - n  •  rpqh  vu  +  |q0 12  v2  +  h 


(6.5) 


We  let 


A(v,h)  =  |q0|2  v2  +h2 ,  B (v,h)  =  2fq0  •  (p0  +  rpqv)  - h  ■  r pqh 


v,  C(v)  =  p0  +  rpqv 


(6.6) 
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We  assume  that  C  is  different  from  zero  and  apply  the  Euler  transformation  ([13],  p.  80) 


R  —  Vc  u  +  Bu  +  A  —  t  —  Cu  . 


(6.7) 


We  square  both  sides  and  solve  for  u 
f-A 

u  = - . 

B  +  2  Ct 

From  this 


(6.8) 


du 


2t(B  +  2Ct)-2C(t2 
(B  +  2  Ctf 


Ct  +  Bt  +  CA 
2 - —dt 

(B  +  20) 


and 


du  1  Ct"  +  Bt  +  CA  2  Ct  +  Bt  +  CA  2  dt 

—  = - 2 - dt  = - , - dt  = - 

R  t-Cu  ( 5  +  20)  f  C  t  -A  ( B  +  2Ct)  B  +  2Ct 

B  +  2Ct 


Applying  the  transformation  to  (6.3),  we  get  that 


C2(l-p0f+B(\-p0)+A+C{l-p0) 

, _ j _ 

^C2(l+p0)2-B(l+p0)+A-C(l+p0) 


2C 

^C2(1-p11)2+B(1-/,„)  +  A+C(1-p„) 

+  B 

2C 

\jc2  (l  +  P0)  -  B(l  + p0)  + A-C(l  + p0) 

+  B 

2dT  =-ln|2a  +  R| 
B  +  2Ct  C  1 


^/c2(l-p0)2+fl(l-p0)+A+C(l-p0 ) 
^C2(l+p0)2-B(l+p0)+A-C(l+p0) 


(6.9) 


(6.10) 


(6.11) 


In  place  of  the  double  integral  in  (6.1),  we  now  have  a  single  integral  with  a  more  complicated 
integrand 


2C -  p0,v)  +  C -  pQ)]  +  B (v,h) 

dv 

2C(v)[r(-(1  +  jp0)a)-C(v)(1  +  p0) 

+  B(v,h ) 

C(v) 

(6.12) 


where  A,  B,  and  C  are  given  by  (6.6). 

We  note  that  the  denominator  of  the  fraction  in  the  algorithm’s  argument  is  equal  to  zero 
when  h  and  v  are  both  zero.  The  numerator,  however,  is  a  positive  number  for  these  values,  as  we 
will  show  at  the  end  of  this  section.  We  thus  have  a  logarithmic  singularity  at  the  origin  when  h 
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is  equal  to  zero  and  the  integrand  goes  to  infinity  at  this  point.  We  turn  to  the  second  BQ  of  Sec¬ 
tion  3  and  use  it  to  test  how  well  Mathematica’s  numerical  subroutines  can  handle  the  integral  in 
(6.12).  The  value  of  po  and  qo  is  as  in  (4.16).  From  (B.2),  and  omitting  the  constant  vector,  we 
find  that  the  IP  is 


r(u,v) 


(  n 

3,0,- 

f  O 
0,2,- 

u  + 

l  4 

l  4 

(0, 0,1)  MV 


(6.13) 


from  which  we  find  that  the  unit  normal  at  the  origin  is 


r 


n  = 


3,0,- 

V  4y 


0,2, 


r 


1 


3,0, 

V  4y 


0,2, 


-2Jr-3y  +  24z 
y/589 


(6.14) 


From  (6.6)  and  (6.13), 


A  65  2  7  2 

A  =  — v  +  h 

16 


B  = 


r  1  96 

v  + - i=h 

V  4  v589  j 


v  „  L 

(  n 

—  ,  C  =  .19  + 

V  H — 

2  V 

l  4) 

(6.15) 


We  substitute  these  results  in  the  integrand  of  (6.12).  In  Figure  6.1,  we  display  the  graph  of 
the  integrand  as  a  function  of  v  over  the  range  of  integration,  [-1.25,  0.75],  for  three  values  of  h. 
We  see  that,  as  h  tends  to  zero,  we  do  get  the  logarithmic  singularity  mentioned  above.  In  terms 
of  numerical  integration,  this  tells  us  that  we  have  a  singularity-type  behavior  at  the  origin  not 
only  when  h  =  0  but  also  when  h  is  a  small,  non-zero  number. 


8 


6 


Figure  6.1.  The  integrand  of  (6.12)  for  the  second  BQ  ( i.e .,  with  (6.15)  substituted  in  it)  as  a 
function  of  v  over  the  range  of  integration.  From  left  to  right:  h  =  0.1,0.01,0.00.  Though  it  is 
not  clear  from  this  picture,  the  function  is  singular  at  the  origin  when  h  =  0. 

We  have  used  two  numerical  integration  routines  in  Mathematica  [11]  to  compute  (6.12)  as 
a  function  of  h.  The  first  method  is  the  GKQ  ([14],  p.  106)  while  the  second  method  is  the  DEQ 
[15],  ([14],  p.  214)  which  is  also  known  as  the  TSQ.  We  exhibit  typical  results  in  Table  6.1.1  and 
6.1.2.  For  each  of  these  methods,  Mathematica  allows  the  specification  of  singular  and  nearly 
singular  points.  In  the  present  case,  we  have  one  such  point.  It  is  the  origin  when  h  is  equal  to  ze¬ 
ro  or  close  to  it.  The  designation  “with”  means  that  we  have  specified  this  point  in  Mathematica, 
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while  the  designation  “w/out”  means  that  we  have  not.  We  set  the  maximum  number  of  recursive 
subdivisions  of  the  integration  interval  to  12  (MaxRecursion  — >  12).  We  also  record  the  execu¬ 
tion  time  in  CPU  sec  and  the  maximum  number  of  SD  we  were  able  to  get  from  each  method.  It 
is  clear  that  specifying  the  (near)  singularity  reduces  execution  time  by  about  an  order  of  magni¬ 
tude.  It  also  provides  more  SD.  We  also  see  that  GKQ  is  faster  than  the  double  exponential 
method  but  the  latter  provides  more  SD.  Whether  the  extra  SD  are  worth  the  time  expense  de¬ 
pends  on  the  application.  Finally,  in  Figure  6.2,  we  display  the  graph  of  the  integral  in  (6.12)  as  a 
function  of  h.  Although  the  integral  appears  to  be  an  even  function  of  h,  it  is  not.  This  is  because 
the  BQ  is  not  flat,  the  only  time  this  property  would  hold.  Simply  put,  the  distance  from  hh  to  a 
point  on  the  BQ  is  not  the  same  as  the  distance  from  —nh  to  the  same  point.  From  these  two  ta¬ 
bles,  we  also  note  that  if  single  precision  (SP)  is  desired,  then  these  algorithms  do  provide  that. 

A  note  of  caution:  the  CPU  times  are  given  for  the  purpose  of  comparison  only,  as  for  ex¬ 
ample  specifying  the  singular  point  vs.  not  specifying  it.  The  computer  on  which  we  ran  these 
calculations  is  a  laptop  with  a  first-generation  Intel  i7  CPU. 

Table  6.1.1.  The  value  of  the  integral  in  (6.12)  for  five  values  of  h  calculated  using 
GKQ  with  and  without  specifying  the  singular  point.  The  execution  time  is  in  CPU 
seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each  meth¬ 
od.  Second  BQ.  Singular  or  nearly  singular  point  at  (p0  =  qo  =  0.25). 


h 

GKQ  w/out 

GKQ  with 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.6488500450514 

0.437 

14 

2.6488500450514 

0.031 

14 

1.00E-04 

2.74952 

0.109 

6 

2.74951658706 

0.047 

12 

1.00E-08 

2.7496 

0.109 

5 

2.7496197 

0.016 

8 

1.00E-12 

2.7496 

0.141 

5 

2.7496196 

0.015 

8 

0.00E+00 

2.7496 

0.109 

5 

2.7496196 

0.031 

8 

Table  6.1.2.  The  value  of  the  integral  in  (6.12)  for  five  values  of  h  calculated  using 
DEQ  with  and  without  specifying  the  singular  point.  The  execution  time  is  in  CPU 
seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each  meth¬ 
od.  Second  BQ.  Singular  or  nearly  singular  point  at  (po  =  qo  =  0.25). 


h 

DEQ  w/out 

DEQ  with 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.6488500450514 

0.53 

14 

2.6488500450514 

0.032 

14 

1.00E-04 

2.7495 

5.35 

5 

2.7495165870543 

0.358 

14 

1.00E-08 

2.7496 

5.289 

5 

2.749619652 

1.201 

10 

1.00E-12 

2.7496 

5.351 

5 

2.74961965 

0.031 

9 

0.00E+00 

2.7496 

5.024 

5 

2.74961965 

0.031 

9 
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Figure  6.2:  Integral  in  (6.12)  as  a  function  of  h. 
We  turn  now  to  the  third  BQ  of  Section  3.  From  (B.18),  we  have  that 
r(«,v)  =  (3,0,2.5)z/  +  (0,2, 2.5);’  +  (0,0,10)wv 
so  that 


10 


v  f 


10 


3,0,—  x  0,2,— 


n  = 


v 


J  V 


(  10") 
3,0,— 
l  4  J 


c 


0,2,  (2 

v  4  j 


-20x-30y +  24£ 

Vl876 


and 


41  , 

A  =  —  v2+/r,  5  =  5 

4 


in  10  96 

10v  + - ; - 

4  V1876 


■,  C  =  ^9  +  (l0v  +  2.5)2 


(6.16) 


(6.17) 


(6.18) 


We  substitute  these  results  in  the  integrand  of  (6.12).  In  Figure  6.3,  we  display  the  graph  of  the 
integrand  as  a  function  of  v  over  the  range  of  integration,  [-1.25,  0.75],  for  three  values  of  h.  The 
behavior  is  the  same  as  for  the  second  BQ.  We  also  compute  the  integral  in  (6.12)  and  form  Ta¬ 
bles  6.2,  tables  analogous  to  Tables  6.1.  The  conclusions  we  draw  from  the  new  tables  are  the 
same  as  those  for  Tables  6.1.  If  we  need  more  than  SP,  a  straight-forward  numerical  evaluation 
of  the  integral  in  (6.12)  is  not  sufficient.  In  the  next  section,  we  proceed  to  re-write  this  integral 
so  as  to  make  it  computable  to  machine  (double)  precision. 
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Figure  6.3.  The  integrand  of  (6.12)  for  the  third  BQ  ( i.e with  (6.18)  substituted  in  it)  as  a  func¬ 
tion  of  v  over  the  range  of  integration.  From  left  to  right:  h  =  0.1,  0.01, 0.00 .  Though  it  is  not 
clear  from  this  picture,  the  function  is  singular  at  the  origin  when  h  =  0. 

Table  6.2.1.  The  value  of  the  integral  in  (6.12)  for  five  values  of  h  calculated  using 
GKQ  with  and  without  specifying  the  singular  point.  The  execution  time  is  in  CPU 
seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each  meth¬ 
od.  Third  BQ.  Singular  or  nearly  singular  point  at  (p0  =  </o  =  0.25). 


h 

GKQ  w/out 

GKQ  with 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

1.8704724171441 

0.453 

14 

1.8704724171441 

0.016 

14 

1.00E-04 

1.92438 

0.109 

6 

1.92437645094 

0.031 

12 

1.00E-08 

1.9244 

0.125 

5 

1.9244343 

0.031 

8 

1.00E-12 

1.9244 

0.156 

5 

1.9244343 

0.015 

8 

0.00E+00 

1.9244 

0.094 

5 

1.9244343 

0.032 

8 

Table  6.2.2.  The  value  of  the  integral  in  (6.12)  for  five  values  of  h  calculated  using 
DEQ  with  and  without  specifying  the  singular  point.  The  execution  time  is  in  CPU 
seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each  meth¬ 
od.  Third  BQ.  Singular  or  nearly  singular  point  at  ( po  =  qo  =  0.25). 


h 

DEQ  w/out 

DEQ  wit 

l 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

1.8704724171441 

0.655 

14 

1.8704724171441 

0.015 

14 

1.00E-04 

1.9244 

5.445 

5 

1.92437645094 

0.187 

13 

1.00E-08 

1.9244 

5.429 

5 

1.924434344 

0.640 

10 

1.00E-12 

1.9244 

5.507 

5 

1.92443434 

0.062 

9 

0.00E+00 

1.9244 

5.116 

5 

1.92443434 

0.047 

9 
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Before  proceeding  further,  we  point  out  that  we  can  also  reverse  the  order  of  integration  in 
(6.1)  and  integrate  first  with  respect  to  v.  In  this  case,  we  get 


2  F(u) 

(m)”  (1  -  ~q0  f  +E  (u,h)(\ -q0)  +  D  ( u,h )  +  F  (u)  (l  -  q0  ) 

+  E ( u,h ) 

clu 

2  F(u) 

\Jf (u)~  (1  +  qQ  f  - E (u,h)(  1  +  q0)  +  D (u,h)  - F (u )(l  +  q0 ) 

+  E ( u,h ) 

F(u) 

(6.19) 


where 


D(«,/?)  =  |p0  u2+h\  E(u,h)-2 


Vo{%+rPqu)-hfl-rPq 


F(u)=  q0+V  •  (6-20) 


Using  (6.19)  is  advantageous  when  the  power  m  in  the  last  integral  in  (5.3)  is  smaller  than  the 
power  /. 

Before  leaving  this  section,  we  remark  that  the  numerator  in  (6.12)  becomes  zero  only  when 
the  singular  point  is  an  edge  point  of  the  BQ.  The  proof  is  as  follows.  Setting  the  numerator 
equal  to  zero,  we  have  that 

^l-p,)2 +«(!-/>„)+ A  =  -^-C(l-p0).  (6.21) 

Ifpo  =  1,  an  edge  point,  then  we  have  that 

41  =  -—.  (6.22) 
2C 


From  (6.6),  this  can  be  satisfied  if,  and  only  if,  v  =  h  =  0. 
If  po  <1,  then,  squaring  both  sides  of  (6.21),  we  get 


A 


V2C, 


From  (6.21)  we  have  that 


B 

- h 

2C 


C(l-p0)<0. 


(6.23) 


(6.24) 


Substituting  (6.23)  in  this,  we  get 

>/A  +  C(l-p0)<0  (6.25) 

which,  by  (6.6),  becomes 
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yj\q0\2v2+h2  +|p0+rMv|(l-p0)<0.  (6.26) 

This  expression  can  never  be  negative.  For  it  to  be  zero,  both  terms  must  be  equal  to  zero.  The 
first  term  is  equal  to  zero  when  both  v  and  h  are  equal  to  zero;  the  second,  however,  is  not.  The 
second  term  becomes  zero  when 


-1 

1 1  —  _ 

TS 

o 

1+ 

[Po‘: 

i2  i  i2 

1  -|Po| 

2 

V  — 

2 

(6.27) 


The  value  of  v  here  is  real  only  when  the  two  vectors  are  collinear;  in  this  case, 


which  makes  the  second  term  in  (6.26)  equal  to  zero  but  not  the  first. 


(6.28) 
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SECTION  7:  EVALUATION  OF  THE  INTEGRAL  IN  (6.12) 

In  this  section  we  develop  a  procedure  for  evaluating  the  integral  in  (6.12)  to  machine  (dou¬ 
ble)  precision. 

As  we  showed  at  the  end  of  Section  6,  the  numerator  of  the  logarithmic  argument  in  (6.12) 
is  never  equal  to  zero  for  interior  singular  points  and,  hence,  it  presents  no  computational  diffi¬ 
culties.  For  this  reason,  we  re-write  (6.12)  in  the  form 


l“9o 


I=  J  ln|2C(v)[/?(l-/?0,v)  +  C(v)(l-/?0)]  +  fl(v,/i)| 


-l-?o 

l-9o 


dv 

cW 


-  J  In  2C(v)[j?(-(l+p„),v)  —  C(v)(l  +  p0)]+B(v,ft) 


dv 


-l-9o 


C(v) 


=  /, 


(7.1) 


where  R(u,  v )  is  defined  in  (6.2).  We  can  evaluate  the  first  integral  numerically  using  either  of 
the  two  methods  of  Section  6.  The  term  in  the  square  brackets  in  the  second  integral  may  become 
unstable  when  both  v  and  h  are  close  to  zero  for,  in  such  a  case,  it  represents  the  difference  of 
two  almost  identical  numbers.  In  Tables  7.1.1  and  7.1.2,  we  have  computed  these  terms  and  their 
difference  at  v  =  0.  We  note  that,  as  h  gets  smaller,  the  value  of  R  approaches  that  of  (l+po)C(v) 
(which  is  independent  of  h )  so  that  their  difference  is  correct  to  very  few  digits.  This  happens  not 
only  at  v  =  0  but  in  an  entire  neighborhood  of  the  origin.  The  BQs  we  used  are  the  second  and 
third  ones  of  Section  3  and  the  BQ  point  is  (0.25,  0.25). 

Table  7.1.1.  The  terms  in  square  brackets  in  the  second  integrand  of  (7.1)  for  v  =  0  and 
for  the  second  BQ.  As  h  tends  to  zero,  the  difference  between  the  two  terms  loses  accu¬ 
racy  and  tends  to  zero,  which  is  the  wrong  answer  if  h  is  different  from  zero. 


h 

^(-(1  +  A))’V) 

(l+p0)C(v) 

Difference 

1.0E-01 

3. 7643 2679904389 E+00 

3. 76299830587259 E+00 

1. 3 28493 17129535 E-03 

1.0E-04 

3. 76299830720132 E+00 

3. 76299830587259 E+00 

1.3 2872779445847 E-09 

1.0E-08 

3. 76299830587259 E+00 

3. 76299830587259 E+00 

0.00000000000000E+00 

1.0E-12 

3. 76299830587259 E+00 

3. 76299830587259 E+00 

0.00000000000000E+00 

1.0E-16 

3. 76299830587259 E+00 

3. 76299830587259 E+00 

0.00000000000000E+00 
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Table  7.1.2.  The  terms  in  square  brackets  in  the  second  integrand  of  (7.1)  for  v  =  0  and 
for  the  third  BQ.  As  h  tends  to  zero,  the  difference  between  the  two  terms  loses  accura¬ 
cy  and  tends  to  zero,  which  is  the  wrong  answer  if  h  is  different  from  zero. 


h 

A(-(l+/?0),v) 

( 1  +po)C(v) 

Difference 

1.0E-01 

4. 88243023503665 E+00 

4.88140604744166E+00 

1. 02418759498679 E-03 

1.0E-04 

4. 88 140604846595 E+00 

4.88140604744166E+00 

1. 0242953 1523285 E-09 

1.0E-08 

4.88140604744166E+00 

4.88140604744166E+00 

0.00000000000000E+00 

1.0E-12 

4.88140604744166E+00 

4.88140604744166E+00 

0.00000000000000E+00 

1.0E-16 

4.88140604744166E+00 

4.88140604744166E+00 

0.00000000000000E+00 

We  correct  for  this  as  follows.  We  let 

8  (v,/*)  =  2C(v) [/?(—(!  +  /70),  v)  —  (1  +  p0)C(v)]  +  5(v,/i) 


From  (6.5) 

R(u,v)  =  p  (v) 
and,  from  (6.6), 

A  =  \q0fv2+h2  , 
For  g  we  write 


+  B(y,ti)u  +  A(y,ti) 


B  =  2 


qo-(po+v;)“"-r/^ 


V , 


C  =  |Po+Vl 


g(v,h)  =  2C(v) 


r2H1+p°)’v)-(1+pq)2c2(v) 

^(-(1  +  Fo)a)  +  (1  +  Fo)C(v) 


(7.2) 


(7.3) 


(7.4) 


=  2  C(v) 


-(1  +  p0)B(v,h)  +  A(v,h) 
R{~(1  +  Po)’V)  +  (l  +  Po)c(v) 


B(v,h ) 


2C(v)[_(1  +  7?o)5(A/?)  +  A(v,/?)]  +  fi(v,/7)[R(-(l  +  jp0),v)  +  (l  +  jp0)C(v)] 

R(-(l  +  Po)’v)  +  (l  +  Po)c(v) 

B(v,h)\_R(-(l  +  p0),v)-(l  +  p0)C(v)]  +  2A(v,h)C(v) 

R(~(l+  Po)a)  +  (1  +  P0)C(v) 
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The  denominator  of  this  fraction  is  never  zero  for  an  interior  point.  In  the  numerator,  the  trouble¬ 
some  term  is  multiplied  by  B.  As  we  can  see  from  (7.4),  B  is  of  0(v),  as  v  tends  to  zero,  and 
sends  the  product  of  the  two  terms  to  zero;  in  fact, 


B(v,h)[R(-(\  +  p0),v)-(\  +  p0)C(v)]  +  2A(v,h)C(v)  ^  >2|p0|/?2  (7.6) 

with  the  contribution  coming  from  the  second  term.  From  (7.5)  and  (7.1),  we  have  that 


Wo 

f  1n 

B(v,h)  7?(-(l  +  p0),v)-(l  +  p0)C(v)  +2 A(v,h)C(v) 

dv 

J  111 

Wo 

^(-(1  +  Po)a)  +  (1  +  P0)C(v) 

C(v) 

Going  back  to  (7.1)  again,  we  write 

1 -<?0  dv 

7=  J  ln|2C(v)[A(l  — p0,v)  +  C(v)(l-p0)]  +  5(v,/7)|— — 

— 1— 9o  V  / 


Wo 

f  in 

B(v,h )  R(-(l  +  p0),v)-(l+  p0)C(v)  +2 A(v,h)C(v) 

dv 

J  ln 

■Wo 

R(-(l  +  Po)’v)  +  (l  +  Po)c(v) 

C(v) 

We  display  results  for  this  expression  in  Tables  7.2.1  (second  BQ)  and  7.2.2  (third  BQ).  For 
both  integration  methods,  we  have  specified  the  origin  as  a  singular  point.  We  have  also  speci¬ 
fied  the  maximum  number  of  recursions  to  be  12.  Only  on  one  occasion  ( h  =  0.1)  do  the  two 
methods  differ  and  only  in  the  last  digit  and  by  one  unit.  The  times  quoted  are  the  CPU  sec  for 
the  two  integrations  in  (7.8).  We  have  not  been  able  to  find  information  in  Mathematica  as  to 
what  zero  time  means.  In  any  case,  the  times  are  shown  here  simply  for  comparing  the  two  inte¬ 
gration  methods.  If  we  compare  these  tables  with  the  right-hand  side  (the  “with”  column)  of  Ta¬ 
bles  6.1  and  6.2,  we  notice  two  things.  First,  that  the  results  of  Tables  6.1  agree  with  the  present 
ones  to  the  number  of  SD  specified  there.  Second,  for  approximately  the  same  time -expense,  we 
get  15  SD  with  the  present  approach. 

The  most  interesting  part  about  the  tables  below  is  the  last  line.  It  clearly  indicates  that  both 
integration  routines  can  handle  even  the  case  h  =  0  and  to  15  SD!  As  we  show  in  the  next  sec¬ 
tion,  this  case  can  be  handled  separately  and  that  the  results  we  will  get  there  are  in  agreement 
with  the  ones  here. 
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Table  7.2.1.  The  value  of  the  expression  in  (7.8)  for  five  values  of  h  calculated  us¬ 
ing  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution  time 
is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from 
each  method.  Second  BQ.  Singular  or  nearly  singular  point  at  (p0  =  q()  =  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.64885004505136 

0.032 

15 

2.64885004505135 

0.046 

15 

1.00E-04 

2.74951658705434 

0.046 

15 

2.74951658705434 

0.031 

15 

1.00E-08 

2.74961966479829 

0.032 

15 

2.74961966479829 

0.063 

15 

1.00E-12 

2.74961967510630 

0.032 

15 

2.74961967510630 

0.046 

15 

1.00E-16 

2.74961967510733 

0.016 

15 

2.74961967510733 

0 

15 

0.00E+00 

2.74961967510733 

0.031 

15 

2.74961967510733 

0 

15 

Table  7.2.2.  The  value  of  the  expression  in  (7.8)  for  five  values  of  h  calculat¬ 
ed  using  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  exe¬ 
cution  time  is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD 
obtainable  from  each  method.  Third  BQ.  Singular  or  nearly  singular  point  at 
(po  =  qo  =  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

1.87047241714406 

0 

15 

1.87047241714406 

0.032 

15 

1.00E-04 

1.92437645094211 

0.031 

15 

1.92437645094211 

0.032 

15 

1.00E-08 

1.92443435241883 

0.047 

15 

1.92443435241883 

0.048 

15 

1.00E-12 

1.92443435820941 

0.063 

15 

1.92443435820941 

0.032 

15 

1.00E-16 

1.92443435820999 

0.047 

15 

1.92443435820999 

0.016 

15 

0.00E+00 

1.92443435820999 

0.030 

15 

1.92443435820999 

0 

15 

The  two  integrals  in  (7.8)  can  be  combined  into  a  single  one 


-% 

f  ln 

{2C(v)[P(l-p0,v)  +  C(v)(l-p0)]  +  B(v,/j)}{i?(-(l  +  p0),v)  +  (l  +  p0)C(v)j 

dv 

J  m 

l-?o 

B(v,h)[R(-(l  +  p0),v)-(l  +  p0)C(v)\  +  2A(v,h)C(v) 

C(v) 

(7.9) 


We  present  results  for  this  integral  in  Tables  7.3.1  and  7.3.2.  The  entries  in  bold  red  indicate 
times  lower  than  those  required  by  (7.8).  From  these  we  conclude  that  GKQ  works  faster  with 
(7.9)  while  for  DEQ  there  is  not  much  time  difference  between  (7.8)  and  (7.9). 
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From  Figs.  6.1  and  6.3,  we  see  that  the  integrand  of  (7.9)  is  not  a  simple  one  to  evaluate 
numerically;  nevertheless,  through  a  simple  manipulation,  we  were  able  to  compute  it  to  machine 
precision  (15  SD). 

Table  7.3.1.  The  value  of  the  integral  in  (7.9)  for  five  values  of  h  calculated  using 
GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution  time  is  in 
CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each 
method.  Second  BQ.  Singular  or  nearly  singular  point  at  (p0  =  qo  =  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.64885004505136 

0.015 

15 

2.64885004505135 

0.031 

15 

1.00E-04 

2.74951658705434 

0.031 

15 

2.74951658705434 

0.015 

15 

1.00E-08 

2.74961966479829 

0.078 

15 

2.74961966479829 

0.062 

15 

1.00E-12 

2.74961967510630 

0.062 

15 

2.74961967510630 

0.032 

15 

1.00E-16 

2.74961967510733 

0.016 

15 

2.74961967510733 

0.015 

15 

0.00E+00 

2.74961967510733 

0.031 

15 

2.74961967510733 

0.016 

15 

Table  7.3.2.  The  value  of  the  integral  in  (7.9)  for  five  values  of  h  calculated  using 
GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution  time  is  in 
CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each 
method.  Third  BQ.  Singular  or  nearly  singular  point  at  (p0  =  q0  =  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

1.87047241714406 

0.015 

15 

1.87047241714406 

0.016 

15 

1.00E-04 

1.92437645094211 

0.015 

15 

1.92437645094211 

0.032 

15 

1.00E-08 

1.92443435241883 

0.031 

15 

1.92443435241883 

0.031 

15 

1.00E-12 

1.92443435820941 

0.016 

15 

1.92443435820941 

0.031 

15 

1.00E-16 

1.92443435820999 

0.016 

15 

1.92443435820999 

0.015 

15 

0.00E+00 

1.92443435820999 

0.016 

15 

1.92443435820999 

0 

15 

In  the  tables  above,  we  notice  that  the  result  for  h  =  0  is  the  same  as  for  h  =  10'16.  We  focus 
on  this  observation  in  Tables  7.4.  There,  we  vary  h  from  10'  to  10'“  and  we  note  that,  for  both 
BQs  and  both  methods,  we  get  the  same  result  as  for  h  =  0  for  values  of  h  equal  or  greater  than 
10‘15.  Considering  (7.6),  and  as  a  precaution,  whenever  h  is  specified  as  zero,  we  may  add  some 
more  stability  to  the  computation  by  giving  it  a  small  value,  such  as  10’“  . 
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Table  7.4.1.  The  value  of  the  integral  in  (7.9)  for  small  values  of  h  calculated  using 
GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution  time  is  in 
CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each 
method.  Second  BQ.  Singular  or  nearly  singular  point  at  (0.25,  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-13 

2.74961967510723 

0.015 

15 

2.74961967510723 

0.016 

15 

1.00E-14 

2.74961967510732 

0.062 

15 

2.74961967510732 

0.016 

15 

1.00E-15 

2.74961967510733 

0.031 

15 

2.74961967510733 

0.015 

15 

1.00E-16 

2.74961967510733 

0.015 

15 

2.74961967510733 

0.016 

15 

1.00E-17 

2.74961967510733 

0.016 

15 

2.74961967510733 

0.015 

15 

1.00E-18 

2.74961967510733 

0.015 

15 

2.74961967510733 

0.000 

15 

1.00E-19 

2.74961967510733 

0.031 

15 

2.74961967510733 

0.000 

15 

1.00E-20 

2.74961967510733 

0.031 

15 

2.74961967510733 

0.000 

15 

0.00E+00 

2.74961967510733 

0.031 

15 

2.74961967510733 

0.000 

15 

Table  7.4.2.  The  value  of  the  integral  in  (7.9)  for  small  values  of  h  calculated  using 
GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution  time  is  in 
CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable  from  each 
method.  Third  BQ.  Singular  or  nearly  singular  point  at  (0.25,  0.25). 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-13 

1.92443435820993 

0.015 

15 

1.92443435820993 

0.016 

15 

1.00E-14 

1.92443435820998 

0.015 

15 

1.92443435820998 

0.016 

15 

1.00E-15 

1.92443435820999 

0.015 

15 

1.92443435820999 

0.000 

15 

1.00E-16 

1.92443435820999 

0.000 

15 

1.92443435820999 

0.015 

15 

1.00E-17 

1.92443435820999 

0.015 

15 

1.92443435820999 

0.000 

15 

1.00E-18 

1.92443435820999 

0.016 

15 

1.92443435820999 

0.015 

15 

1.00E-19 

1.92443435820999 

0.000 

15 

1.92443435820999 

0.016 

15 

1.00E-20 

1.92443435820999 

0.015 

15 

1.92443435820999 

0.000 

15 

0.00E+00 

1.92443435820999 

0.016 

15 

1.92443435820999 

0.015 

15 

In  the  next  section,  we  present  more  examples  and  another  way  of  evaluating  (7.9)  when 
h  =  0. 
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SECTION  8:  ADDITIONAL  TEST  CASES 

In  the  previous  sections,  we  developed  a  method  for  computing  the  inner  integral  of  the  IM 
element  of  EFIE.  We  tested  the  method  using  the  BQs  of  Section  3  but  we  employed  only  one 
point  in  the  parametric  representation  of  each  BQ,  namely,  the  point  po  =  qo  =  0.25.  In  this  sec¬ 
tion,  we  expand  the  testing  using  the  third  BQ  of  Section  3.  As  we  pointed  out  in  that  section, 
this  extreme  type  of  a  BQ  should  be  avoided  in  practice.  It  makes,  however,  for  a  good  test  case 
since,  if  our  method  can  provide  good  results  for  this  BQ,  it  will,  in  all  likelihood,  provide  good 
results  for  more  commonly  used  BQs,  such  as  the  second  one  in  Section  3. 

The  first  case  we  consider  is  the  case  po  =  qo  =  0.  At  this  point  the  BQ  has  a  saddle  type  of 
behavior.  From  (3.9),  the  equation  for  the  third  BQ  is 

r{p,q)  =  (3,0,0) p  +  (0,2,0)$  +  (0,0,10) pq  ,  |^]<1,  \q\  <  1 
from  which  we  get  that 

‘3l^  =  (3,0,0)  +  (0,0,10)</,  /1^  =  (0,2,0)  +  (0,0,10)p 
op  oq 


(8.1) 


(8.2) 


and 


(3, 0, l(k/)x(0, 2,1  Op)  _  (3,0,10g)x(0,2,10p)  _  -20qx-30py  +  6z 
"■  |(3,0,10?)x(0,2,10/>)|  ”  1(3,0, 10?)x(0,2,10/>)|  “  ^40092  +  900p2  +36  '  '  ' 

We  use  this  information  in  the  integral  under  consideration,  namely  (7.9).  We  display  the 
graph  of  the  integrand  for  h  =  0.1  in  Figure  8.1.  In  Table  8.1,  we  display  the  results  of  the  inte¬ 
gration  of  (7.9)  using  both  GKQ  and  DEQ,  and  for  the  same  values  of  h  that  we  have  used  in  the 
previous  sections.  The  conditions  for  the  two  quadratures  are  the  same  as  before,  with  the  num¬ 
ber  of  recursions  limited  to  12.  We  have  also  specified  the  origin  as  a  singular  point.  We  obtain 
agreement  to  15  SD,  except  in  the  last  two  cases  where  there  is  a  disagreement  of  one  unit  in  the 
last  digit. 
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Figure  8.1.  The  graph  of  the  integrand  in  (7.9)  for  the  third  BQ  of  Section  3,  with 
Po  =  Qo  =  0  and  h  =  0.1.  As  h  gets  smaller,  the  peak  moves  upward.  When  h  =  0, 
we  have  a  logarithmic  singularity  at  the  origin. 

Table  8.1.  The  numerical  integration  of  (7.9)  for  the  third  BQ  of  Section  3  and  with 
Po  =  qo  =  0.  Both  methods  agree  except  for  the  numbers  in  red,  which  disagree  in 
the  last  digit  by  one  unit. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.16084390813449 

0.047 

15 

2.16084390813449 

0.047 

15 

1.00E-04 

2.25593244874981 

0.015 

15 

2.25593244874981 

0.047 

15 

1.00E-08 

2.25603714690303 

0.031 

15 

2.25603714690303 

0.031 

15 

1.00E-12 

2.25603715737396 

0.031 

15 

2.25603715737396 

0.031 

15 

1.00E-16 

2.25603715737501 

0.032 

14 

2.25603715737500 

0.015 

14 

0.00E+00 

2.25603715737501 

0.016 

14 

2.25603715737500 

0.016 

14 

We  next  give  the  singular  point  the  value  po  =  -1,  qo  =  1.  The  position  vector  from  (8.1)  is  (- 
3,  2,  -10)  and  corresponds  to  one  of  the  two  lower  comers  of  the  BQ.  In  Figure  8.2,  we  show  the 
graph  of  the  integrand  of  (7.9)  for  h  =  0.1.  The  same  comments  as  above  hold  with  the  exception 
that  the  logarithmic  singularity  is  now  an  endpoint  of  the  integration  interval.  We  have  calculated 
the  integral  in  (7.9)  and  we  display  the  results  in  Table  8.2.  We  see  that  this  time  the  two  quadra¬ 
tures  differ  for  large  rather  than  small  values  of  h.  Since  the  difference  is  only  in  the  last  digit 
and  only  by  a  unit,  we  do  not  think  it  is  a  matter  to  dwell  on. 
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Figure  8.2.  The  graph  of  the  integrand  in  (7.9)  for  the  third  BQ  of  Section  3,  with 
Po  =  -1,  qo  =  1  and  h  =  0.1.  As  h  gets  smaller,  the  peak  moves  upward.  When  h  = 
0,  we  have  a  logarithmic  singularity  at  the  origin. 

Table  8.2.  The  numerical  integration  of  (7.9)  for  the  third  BQ  of  Section  3  and  with 
Po  =  -1,  qo  =  1.  Both  methods  agree  except  for  the  numbers  in  red,  which  disagree 
in  the  last  digit  by  one  unit. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

0.421015601245452 

0.031 

15 

0.421015601245452 

0.032 

15 

1.00E-04 

0.422730809605921 

0.031 

14 

0.422730809605920 

0.031 

14 

1.00E-08 

0.422732511450675 

0.031 

14 

0.422732511450674 

0.016 

14 

1.00E-12 

0.422732511620856 

0.032 

15 

0.422732511620856 

0.031 

15 

1.00E-16 

0.422732511620873 

0.015 

15 

0.422732511620873 

0.016 

15 

0.00E+00 

0.422732511620873 

0.016 

15 

0.422732511620873 

0.031 

15 

We  chose  the  next  singular  point  to  be  at  p0  =  <?o  =  1.  The  position  vector  from  (8.1)  is  (3,  2, 
10)  and  corresponds  to  one  of  the  two  upper  comers  of  the  BQ.  In  Figure  8.2,  we  show  the  graph 
of  the  integrand  of  (7.9)  for  h  =  0.1.  We  have  calculated  the  integral  in  (7.9)  and  we  display  the 
results  in  Table  8.2.  We  see  that  this  time  the  two  quadratures  differ  in  all  cases  but  one,  but  only 
in  the  last  digit  and  not  more  than  by  a  unit.  From  the  last  three  cases,  we  may  conclude  that,  in 
terms  of  the  Mathematica  quadratures,  our  method  provides  at  least  14  SD. 
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Figure  8.3.  The  graph  of  the  integrand  in  (7.9)  for  the  third  BQ  of  Section  3,  with 
Po  =  Qo  =  1  and  h  =  0.1.  As  h  gets  smaller,  the  peak  moves  upward.  When  h  =  0, 
we  have  a  logarithmic  singularity  at  the  origin. 

Table  8.3.  The  numerical  integration  of  (7.9)  for  the  third  BQ  of  Section  3  and  with 
Po  =  Qo  -  1.  The  numbers  in  red  indicate  a  disagreement  in  the  last  digit  of  one  unit. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

0.422500945889474 

0.016 

15 

0.422500945889474 

0.015 

15 

1.00E-04 

0.422732296117918 

0.016 

14 

0.422732296117917 

0.015 

14 

1.00E-08 

0.422732511599326 

0.015 

14 

0.422732511599325 

0.016 

14 

1.00E-12 

0.422732511620871 

0.062 

14 

0.422732511620870 

0.015 

14 

1.00E-16 

0.422732511620873 

0.016 

14 

0.422732511620872 

0.015 

14 

0.00E+00 

0.422732511620873 

0.031 

14 

0.422732511620872 

0 

14 

We  also  repeat  the  calculations  as  in  Tables  7.4  for  the  three  points  of  this  section.  The  con¬ 
clusions  are  the  same  as  in  Section  7:  giving  h  a  small  value  when  it  is  supposed  to  be  zero  may 
stabilize  the  algorithm  without  altering  the  end  result. 
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Table  8.4.1.  The  value  of  the  integral  in  (7.9)  for  the  third  BQ  and  for  small  values 
of  h  calculated  using  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The 
execution  time  is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD  ob¬ 
tainable  from  each  method.  Singular  or  nearly  singular  point  at  po  =  0.0,  qo  =  0.0. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-13 

2.25603715737490 

0.047 

15 

2.25603715737490 

0.016 

15 

1.00E-14 

2.25603715737500 

0.031 

14 

2.25603715737499 

0.016 

14 

1.00E-15 

2.25603715737501 

0.000 

14 

2.25603715737500 

0.016 

14 

1.00E-16 

2.25603715737501 

0.032 

14 

2.25603715737500 

0.031 

14 

1.00E-17 

2.25603715737501 

0.031 

14 

2.25603715737501 

0.016 

14 

1.00E-18 

2.25603715737501 

0.016 

14 

2.25603715737500 

0.015 

14 

1.00E-19 

2.25603715737501 

0.031 

14 

2.25603715737501 

0.015 

14 

1.00E-20 

2.25603715737501 

0.031 

14 

2.25603715737500 

0.016 

14 

0.00E+00 

2.25603715737501 

0.031 

14 

2.25603715737501 

0.016 

14 

Table  8.4.2.  The  value  of  the  integral  in  (7.9)  for  the  third  BQ  and  for  small  values 
of  h  calculated  using  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The 
execution  time  is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD  ob¬ 
tainable  from  each  method.  Singular  or  nearly  singular  point  at  po  =  -1.0,  qo  =  1.0. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-13 

0.422732511620871 

0.015 

15 

0.422732511620871 

0.032 

15 

1.00E-14 

0.422732511620873 

0.016 

15 

0.422732511620873 

0.031 

15 

1.00E-15 

0.422732511620873 

0.032 

15 

0.422732511620873 

0.015 

15 

1.00E-16 

0.422732511620873 

0.016 

15 

0.422732511620873 

0.031 

15 

1.00E-17 

0.422732511620873 

0 

15 

0.422732511620873 

0.032 

15 

1.00E-18 

0.422732511620873 

0.031 

15 

0.422732511620873 

0.015 

15 

1.00E-19 

0.422732511620873 

0.031 

15 

0.422732511620873 

0.016 

15 

1.00E-20 

0.422732511620873 

0.015 

15 

0.422732511620873 

0.000 

15 

0.00E+00 

0.422732511620873 

0.016 

15 

0.422732511620873 

0.031 

15 
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Table  8.4.3.  The  value  of  the  integral  in  (7.9)  for  the  third  BQ  and  for  small  values 
of  h  calculated  using  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The 
execution  time  is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD  ob¬ 
tainable  from  each  method.  Singular  or  nearly  singular  point  at  po  =  1.0,  c/u  =  1.0. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-13 

0.422732511620873 

0.016 

14 

0.422732511620872 

0.015 

14 

1.00E-14 

0.422732511620873 

0.046 

14 

0.422732511620872 

0.016 

14 

1.00E-15 

0.422732511620873 

0.016 

14 

0.422732511620872 

0.016 

14 

1.00E-16 

0.422732511620873 

0.031 

14 

0.422732511620872 

0.016 

14 

1.00E-17 

0.422732511620873 

0.000 

14 

0.422732511620872 

0.031 

14 

1.00E-18 

0.422732511620873 

0.031 

14 

0.422732511620872 

0.000 

14 

1.00E-19 

0.422732511620873 

0.000 

14 

0.422732511620872 

0.016 

14 

1.00E-20 

0.422732511620873 

0.016 

14 

0.422732511620872 

0.015 

14 

0.00E+00 

0.422732511620873 

0.031 

14 

0.422732511620872 

0.000 

14 

We  next  present  another  way  of  computing  (7.9)  when  h  is  equal  to  zero.  We  repeat  here 
(7.9)  for  convenience 


!-?o 


/=  |  in 


-!-?o 


[2C(v)[R(l- p0,v)  +  C(v)(l- p0)]  +  B(v,h)}{R(-(l  +  p0),v)  +  (\  +  pQ)C(v)} 


B(v,/i)[i?(-(l  +  p0),v)-(l  +  /?0)C(v)]  +  2A(v,/r)C(v) 


dv 


where,  from  (7.3), 

R (u, v)  =  ^C2  (v)m2  +  B(v,h)u  +  A(v,h) 
while,  from  (7.4), 

^(v,0)  =  |q0|V,  £(v,0)  =  2q0-(p0+rp9v)v,  C(v)  =  |p0  +rpqv 


C(v) 

(8.4) 

(8.5) 


(8.6) 


We  consider  the  case  of  an  interior  singular  point.  When  h  =  0,  the  numerator  of  the  loga¬ 
rithm’s  argument  does  not  become  zero,  nor  does  it  become  so  at  any  other  value  of  the  integra¬ 
tion  variable.  For  the  denominator,  we  write 


(v)  =  B(v,  °)  [/?(—(!  +  />o)>v)  —  (1+7,o)C(v)]  +  2A(v,  °)C(v) 
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=  B(v,  0) 


Ri(-(l  +  Po),v)-(l  +  Po)2C2(v) 


+  2A(v,0)C(v) 


^(-(1  +  Po)a)  +  (1  +  P0)C(v) 

5 (v, 0)  /?2  (- (1  +  P0 ) , v)  -  (1  +  p0 )2  C2  (v)  +  2A(v, 0) C (v) [R (- (1  +  pQ ) , v)  +  (l  +  p0 ) C (v) 


^(-(1  +  Po)a)  +  (1  +  R0)c(v) 


(8.7) 


We  evaluate 


^2(-(1  +  ^o)’v)-(l  +  P0)2  c~  {v)  =  (l  +  Po)2  C2(v)  +  B(v,0)u  +  A(v,0)-(l  +  poy  C 2  (v) 

=  -B  (v,  0)  (1  +  p0 )  +  A(v,  0)  =  -2(l  +  p0)  q0  •  (p0  +  rpcv)  v  +  |q0 12  v2  (8.8) 


and  substitute  in  (8.7) 

5(v,0)  -2(l  +  p0)q0-(p0+rwv)v  +  |q0|V  +2A(v,0)C(v)[/?(-(l  +  p0),v)  +  (l  +  p0)C(v)] 


s(v)  = 


—  2v2 


B(~(l+  Po)’v)  +  (1  +  Po)C(v) 

qo-(Po+v,)[“2(1+Po)qo-(Po+v')+l(ior^]+|qor<7(v)[R(-(i+p0),v)+(i+p0)c(v)] 


^(-(1  +  Po)’v)  +  (1+Po)C(v) 


(8.9) 


We  set  h  =  0  in  the  argument  of  the  logarithm  in  (8.4)  and  write 


in 


N(v) 


2  v2D(v) 


=  ln 


N(v) 


2  D(v) 


-21n  v 


(8.10) 


where 


7V(v)  =  {2C(v)[/?(1  —  Po>v)  +  C(v)(l  —  P0)]  +  5(v,0)}{/?(- (l  +  jp0),v)  +  (l  +  jp0)C(v)}2 

(8.11) 


and 


D  (v)  =  q„  ’  (P0  +  rPqv)  [-2(1  +  A, )  (p0  +  rpqv)  •  q0  +  |q 
+|q0|2c(v)[/?(-(i+jp0),v)+(i+jp0)c(v) 


(8.12) 


We  return  to  (8.4)  and  write 
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l-^0  1 

f  1  ln 

7V(v) 

'7°  i 

✓/i  ’  —  1  In 

JV(v) 

rfv-2?  '"H 

Jam 

2v2D(v) 

civ  —  m 

-f-AM 

2D(v) 

Jaw 

-dv-Ix  (0)-2/2  (0). 

(8.13) 

We  note  that  the  argument  of  the  logarithm  in  the  first  integral  on  the  right  takes  on  a  finite  value 
when  v  is  equal  to  zero.  For  the  second  integral,  we  have 


.  ,  '7“  lnlvl  °r  lnlvl  '7’  In (v)  *7°  ln(v)  '7°ln(v) 

We  integrate  the  new  integrals  once  by  parts 

(i+^0)[in(i+4,o)-1]  (Wo^M1  -%)-!] 


i2(o) 


C{-\-q0) 

l+9o  V  D»  (4-1] 
+  1  - 


C(l~q0 ) 


r«  v-Po’rp, 


C3(-v) 


dv+  | 


?o  v[ln(v)-l] 


lrJ  v  +  Po’rp9 


(8.14) 


C3(v) 


dv.  (8.15) 


The  last  two  integrals  have  no  singularities  and  are  stable  for  numerical  evaluation.  If  we  substi¬ 
tute  this  result  in  (8.14),  we  get 


'«>)  =  1 


1 


Lc(v) 


In 


N(v) 


2  D{v) 


dv  —  2  ■ 


(!  +  )  [ln  (!  +  )  -  i]  (i-^ojpn^-^o)-1] 


C(-l-*o) 


-2 


ijo  v[ln(v)-l] 


y-Po'r« 


C3 (-v) 


*7°  v[ln(v)  - 1] 
dv+  J  - 


C{l-q0) 


T,J  v  +  p0-rp 


C3(v) 


dv 


.(8.16) 


All  the  integrals  in  this  expression  are  numerically  stable. 

We  examine  next  the  case  in  which  po  takes  on  its  extreme  values.  When  p0  -  -1 ,  we  still 

have  the  denominator  of  the  logarithmic  function  becoming  zero  but  not  the  numerator.  In  this 
case, 


^2(°’v)U=o  =|q< 


IV 


A2  (2, 


=  4C2(v)  +  25(v,0)  +  A(v,0);  (^0=-l) 


(8.17) 


and 


2  v2Z)(v) 


2  C(v) 

^4C2  (v)  +  25(v,0)  +  A(v,0)  +  2  C(v) 

+  5(v,0) 

2v 

q0-l 

(po+v)+  q0 

C(vf 

(8.19) 
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If  we  take  v  away  from  the  denominator,  then  the  rest  of  the  fraction  attains  a  finite  value  at  v  = 
0.  We  can  then  proceed  as  above. 

When  p0  =  1 ,  we  have  from  the  end  of  Section  6  that  the  numerator  in  the  logarithm's  argu¬ 
ment  becomes  zero.  In  this  case,  we  have  that 


JV(v) 

{c  (v)  |q0 1  +  q0  •  (p0  +  rpqv)  (v, 0)  +  2C  (v) 

f 

2v-D(v)  „ 

q0-(po+v) 

-4(po+v)-q0+|qor 

*  +  q02c(v) 

i  i 

In 

o 

+ 

N> 

i  i 

(8.20) 

with  same  comments  as  above.  The  singularity  is  of  the  same  kind  as  the  one  in  (8.19).  We  recall 
that  we  used  the  last  two  values  of  po  in  the  examples  above  and  all  cases  ran  smoothly. 

We  return  to  (8.16)  and  we  compute  it  for  the  third  BQ  and  the  same  points.  The  results  are 
displayed  in  Table  8.5.  The  results  of  the  first  row  are  in  agreement  with  those  of  Table  7.3.2 
(last  row).  The  results  of  the  second  row  are  in  agreement  with  those  of  Table  8.3.1  (last  row)  but 
took  a  lot  longer.  The  results  of  the  third  row  are  in  agreement  with  those  of  Table  8.3.2  (last 
row).  From  the  CPU  time,  we  see  that  GKQ  struggled  to  get  to  the  result  while  the  DEQ  time  is 
the  same  as  in  Table  8.3.2.  From  the  last  row,  we  see  that  both  quadratures  achieved  15  SD  but 
the  CPU  time  they  required  is  much  higher  than  that  listed  in  Table  8.3.3.  The  conclusion  here  is 
that  we  will  do  as  well  using  (7.9)  instead  of  (8.17). 

Table  8.5.  The  value  of  (8.17)  for  the  third  BQ  and  for  four  values  of  ( po ,  t/o)  calcu¬ 
lated  using  GKQ  and  DEQ  methods  and  specifying  the  singular  point.  The  execution 
time  is  in  CPU  seconds.  We  also  display  the  maximum  number  of  SD  obtainable 
from  each  method. 


(po,  do) 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

(0.25,  0.25) 

1.92443435820999 

0.015 

15 

1.92443435820999 

0.016 

15 

(0.0,  0.0) 

2.25603715737501 

0.046 

15 

2.25603715737500 

0.046 

15 

(-1.0, 1.0) 

0.422732511620873 

0.110 

15 

0.422732511620873 

0.030 

15 

(1.0,  1.0) 

0.422732511620873 

0.093 

15 

0.422732511620873 

0.048 

15 

We  conclude  this  section  with  a  brief  explanation  as  to  why  some  of  the  times  in  Table  8.4 
are  much  higher  than  the  corresponding  times  in  the  previous  tables.  In  Figure  8.4,  we  display 
the  integrand  of  the  middle  integral  in  (8.16)  for  po  =  -1  and  do  =  1 .  The  fact  that  we  are  subtract¬ 
ing  two  almost  equal  areas  may  have  contributed  to  the  slow  convergence  of  GKQ.  The  execu¬ 
tion  time  in  this  case  is  0.094  CPU  sec  vs.  0.015  for  the  DEQ.  The  same  is  true  for  the  case  po  = 
1  and  c/t)  =  1 .  We  display  the  corresponding  integrand  in  Figure  8.5. 
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Figure  8.4:  Integrand  of  the  middle  integral  in  (8.16)  for  po  =  -1  and  qo 


V 


Figure  8.5:  Integrand  of  the  middle  integral  in  (8.16)  for  p0  =  1  and  q0 
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SECTION  9:  PROOF  OF  CFAIM  AND  SENSITIVITY  ANAFYSIS 

In  Section  6,  we  claimed  that  C  in  (6.6)  is  greater  than  zero.  We  show  below  that,  if  (6.6)  is 
equal  to  zero,  we  have  an  irregular  BQ.  We  also  give  additional  examples  of  irregular  BQs  and 
perform  a  sensitivity  study  to  determine  how  precision  is  affected  as  we  pass  from  a  regular  to  an 
irregular  BQ.  The  arguments  that  we  use  here  can  also  be  used  to  show  that  F  in  (6.20)  is  also 
greater  than  zero.  From  (6.6),  we  have  that 

C  =  |p„+rwv|.  (9.1) 

We  assume  first  that  rpq  =  0 .  From  (2.2),  we  see  that  the  resulting  BQ  describes  a  planar 
parallelogram.  For  (9.1)  to  be  zero,  we  must  have  that  p()  =  0 .  From  (2.23)  the  position  vector  to 
the  BQ  is 


r(p,  q )  =  r0  +  (r,  +  r pqq0 ){p~p0)  +  (r  +  rpqPo ){q-q0)  +  rpq  (p  -  pQ  )(q-q0)  (9.2) 

with  ro  defined  in  (2.22).  From  (2.25) 


„  _5r(p0,d0)  „  _dr(p0,qQ ) 

Po  ~  >  4o  ~  ~  \^-J) 

OP0  oqQ 

so  that,  from  the  last  two, 

Po=rP+rPq%-  (9-4) 

We  let  now  p0  and  rpq  be  zero  in  (9.2)  and  we  get 

r(p,q)  =  r0+rg(q-q0)  (9.5) 

which  is  the  equation  of  a  straight  line  in  space.  We  thus  have  an  irregular  BQ.  We  could  have 
surmised  that  immediately  since  p0  =  0  implies  that  (2.7)  is  zero. 

We  assume  next  that  rpq  ^0 .  For  (9.1)  to  be  equal  to  zero,  we  must  have  the  two  vectors 

there  to  be  collinear.  From  this  and  (9.4),  we  conclude  that  rpq  and  rp  are  collinear.  This  implies 
that  the  BQ  is  planar.  We  write 

rpq  =  ccrp  .  (9.6) 

If  we  compare  this  expression  with  (2.9),  we  conclude  that  (3=  0  and  Table  2.1  tells  us  for  what 
values  of  a  we  will  have  a  regular  or  irregular,  planar  BQ.  We  verify  these  values  by  substituting 
(9.6)  in  (9.1) 
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C  =  P0  +  r  v  =  rp  +r pqq0  +rpqv  =  rp  +r ■  q  =  rp  +r paq  =  rp  \l  +  aq\ 


(9.7) 


For  this  to  be  equal  to  zero,  we  must  have  that 

1 


q  =  — 


a 


(9.8) 


If  \a\  <  1 ,  the  absolute  value  of  q  that  makes  (9.7)  be  zero  is  greater  than  one.  Since  -  \<q<\, 
then  (9.7)  is  different  from  zero.  If  \a\  >  1 ,  then  \q\  <  1  and  we  can  have  (9.7)  equal  to  zero.  What 
does  this  imply  in  terms  of  the  nature  of  the  BQ?  If  we  substitute  (9.6)  in  (9.2),  we  get 

r (P, q)  =  r0  +  (r,  +  arpqQ )(p-p0)  +  (rg+  arpPo )(q-q0)  +  arp  ( p  -  p0 )(q~q0) 

=  rn  +  rP  [(!  +  ccq0)(p-  p0  )  + ap0(q  -  q0)  + a(p  -  p0  )  (q  ~  q0  )]  +  rq  (q  -  q0  )  .  (9.9) 

From  this,  we  see  that  we  have  a  surface  on  the  plane  defined  by  the  vectors  rp  and  rq,  as  we 
have  already  mentioned.  We  proceed  to  compute 


dr(p,q) 

dp 


i[(l  +  aqQ)  +  a(q-qQ)], 


dr(p,q) 

dq 


rpa(p-pt))  +  rr 


(9.10) 


From  this 


dr(p,q)  dr(p,q) 

- X - 


dp 


dq 


=  [(1  +  aqQ )  ■ +  a  (q  -  q0 )] (rp  x  ^ rq )  -  (l  +  aq)  (rp  x  rq ) 


(9.11) 


which  is  equal  to  zero  when  q  is  given  by  (9.8).  These  results  are  in  agreement  with  Table  2.1. 
Thus,  the  only  time  (9.1)  can  be  zero  is  when  we  have  an  irregular  BQ. 

We  have  given  examples  of  regular  and  irregular  BQs  in  Section  2.  We  present  additional 
ones  here  by  appropriately  modifying  the  third  BQ  of  Section  3.  The  relevant  vectors  are 

rp  =(3,0,0),  r,  =(0,2,0),  r^=(3«,0,0).  (9.12) 


We  have  thus  modified  the  last  vector  to  be  collinear  with  rp.  We  also  set  po  =  qo  =  0.25.  Since 
the  vectors  in  (9.12)  lie  in  the  xy-plane,  so  does  the  surface  of  the  BQ.  In  Figure  9.1.1,  we  display 
the  BQ's  surface  when  a  =  0.5.  Since  (9.8)  results  in  a  value  of  q  outside  its  range  of  values,  we 
have  a  regular  BQ  in  the  xy-plane.  Moreover,  for  C  in  (9.1)  we  have  that 


C  = 


)0+V 


rp(l  +  aq)\  =  3 


1  + 


q 


> 


(9.13) 
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In  Figure  9.1.2,  we  have  set  a  =  1.0.  The  resulting  value  of  -1  in  (9.8)  is  at  the  lower  end  of  the 
range  of  q  and  results  in  the  QL  being  reduced  to  a  triangle: 

C  =  \rp(l  +  aq)\=3(l  +  q)>0.  (9.14) 

In  Figure  9.1.3,  we  have  set  a  =  2.0.  The  resulting  value  of  q  in  (9.8)  is  now  well  within  the 
range  of  values  of  q  and  we  get  a  very  irregular  BQ,  namely,  two  triangles  with  a  common  ver¬ 
tex. 


-3  <  C  =  3(1  +  2q)  <  9 . 

We  can  obtain  similar  results  for  F  in  (6.20). 


(9.15) 


Figure  9.1.1:  Graph  of  the  BQ  of  (9.12)  with  a=  0.5. 


Figure  9.1.2:  Graph  of  the  BQ  of  (9.12)  with  a=  1.0. 


65 


NAWCADPAX/TR-20 15/214 


2 


0 


x 


-1 


-2 


-5 


0 


5 


Figure  9.1.3:  Graph  of  the  BQ  of  (9.12)  with  a=  2.0. 


In  this  last  case,  the  position  vector  for  the  BQ  is 


r(p,q)  =  xp(\  +  2q)+y2q. 


(9.16) 


From  this  we  see  that  the  line  q  =  -  1/2  ( |/;|  <  1)  maps  to  the  single  point  (0,  -1)  in  the  xy-plane. 

From  the  above  discussion,  it  becomes  clear  that  we  have  to  exercise  caution  when  we  deal 
with  a  flat  BQ.  From  (2.2),  the  position  vector  to  a  point  on  the  BQ  is 


(9.17) 


r(p,  q)  =  r00  +  rpp  +  rqq  +  rpqpq  ,  \p\<l,  \q\  <  1 


with  the  four  vectors  defined  in  (2.3)  -  (2.5).  For  this  BQ  to  be  flat,  we  must  have  the  last  three 
vectors  lying  on  the  same  plane.  One  way  to  express  this  condition  is  that 


(9.18) 


If  this  holds,  the  next  step  is  to  determine  the  scalars  in  (2.9).  This  can  be  done  by  considering 
the  rectangular  components  of  this  equation  and  solving  a  system  of  two  equations  in  two  un¬ 
knowns.  Once  the  unknowns  have  been  determined,  we  can  use  Table  2.1  to  decide  whether  we 
are  dealing  with  a  regular  or  irregular  BQ. 

We  next  ask  the  question  of  how  sensitive  the  numerical  integration  becomes  as  we  move 
from  a  regular  BQ  to  an  irregular  one.  To  this  end,  we  introduce  the  BQ 


rp=  (3,0,0),  r?  =(0,2,0),  rpq  =(0,2j3,0) . 


(9.19) 


For  this  BQ,  rq  and  rpq  are  collinear.  In  place  of  (9.6),  we  now  have 


r 


(9.20) 
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while,  in  place  of  (9.8),  we  have 


In  Figures  9.2.1  to  9.2.3,  we  display  the  graph  of  this  BQ  when  /?  =  0.9,  0.99  and  0.999.  We  see 
that,  as  /?  tends  to  one,  the  BQ  tends  to  look  more  like  a  triangle  than  a  QL. 


y 


Figure  9.2.1:  Graph  of  the  BQ  of  (9.19)  with  / 3  =  0.9  ip'  =  -1.1111). 


y 


Figure  9.2.2:  Graph  of  the  BQ  of  (9.19)  with  /?=  0.99  ip'  =  -1.0101). 
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y 


Figure  9.2.3:  Graph  of  the  BQ  of  (9.19)  with  (3  =  0.999  ip'  =  -1.0010). 

We  now  examine  what  happens  to  the  integral  in  (7.9)  as  /?  tends  to  one.  We  have  computed 
this  integral  for  the  singular  point  at  the  center  of  the  /^/-coordinates  and  also  at  a  comer  point. 
We  have  considered  the  three  values  of  [3  as  well  as  four  values  of  h.  We  display  the  integration 
results  in  Tables  9.1.1  to  9.2.3.  We  see  that  the  two  quadratures  agree  to  at  least  14  SD  and, 
where  there  is  a  disagreement  in  the  last  digit,  it  is  only  by  one  unit.  We  also  note  that  the  cases 
of  h  =  10''  and  h  =  0  are  in  perfect  agreement,  which  further  reinforces  the  use  of  a  small,  non¬ 
zero  h  in  place  of  zero. 

Table  9.1.1:  Value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  qo 
=  0  and  (3  =  0.9. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.74665336952689 

0.016 

14 

2.74665336952688 

0.031 

14 

1.00E-08 

2.84965918230493 

0.031 

15 

2.84965918230493 

0.032 

15 

1.00E-20 

2.84965919277691 

0.031 

14 

2.84965919277690 

0.015 

14 

0.00E+00 

2.84965919277691 

0.031 

14 

2.84965919277690 

0.000 

14 

Table  9.1.2:  Value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  qo  = 
0  and  /3=  0.99. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.74730724784908 

0.000 

15 

2.74730724784908 

0.016 

15 

1.00E-08 

2.85037159343782 

0.032 

14 

2.85037159343781 

0.031 

14 

1.00E-20 

2.85037160390979 

0.016 

15 

2.85037160390979 

0.015 

15 

0.00E+00 

2.85037160390979 

0.015 

15 

2.85037160390979 

0.015 

15 
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Table  9.1.3:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  qo 
=  0  and  /?=  0.999. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.74736479069567 

0.000 

14 

2.74736479069566 

0.016 

14 

1.00E-08 

2.85043525566924 

0.031 

15 

2.85043525566924 

0.031 

15 

1.00E-20 

2.85043526614122 

0.032 

15 

2.85043526614122 

0.015 

15 

0.00E+00 

2.85043526614122 

0.016 

15 

2.85043526614122 

0.016 

15 

Table  9.2.1:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  = 
-1,  qo  =  1,  and  J3  =  0.9. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.56501461685842 

0.015 

15 

2.56501461685842 

0.032 

15 

1.00E-08 

2.84064664644604 

0.031 

15 

2.84064664644604 

0.016 

15 

1.00E-20 

2.84064668163297 

0.016 

15 

2.84064668163297 

0.015 

15 

0.00E+00 

2.84064668163297 

0.000 

15 

2.84064668163297 

0.000 

15 

Table  9.2.2:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  - 
1,  qo  =  1,  and  J3  =0.99. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.99932746262052 

0.016 

14 

2.99932746262051 

0.015 

14 

1.00E-08 

4.26080553788516 

0.015 

15 

4.26080553788516 

0.016 

15 

1.00E-20 

4.26080589691273 

0.000 

15 

4.26080589691273 

0.016 

15 

0.00E+00 

4.26080589691273 

0.016 

15 

4.26080589691273 

0.015 

15 

Table  9.2.3:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  = 
-1,  qo  =  1,  and  J3  =0.999. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

3.02996129711238 

0.000 

15 

3.02996129711238 

0.000 

15 

1.00E-08 

5.69681413172724 

0.000 

14 

5.69681413172723 

0.016 

14 

1.00E-20 

5.69681772889930 

0.016 

14 

5.69681772889929 

0.015 

14 

0.00E+00 

5.69681772889930 

0.015 

14 

5.69681772889929 

0.016 

14 
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We  may  also  ask  the  question  of  whether  (7.9)  is  computable  when  \j3\  >  1 .  Here,  we  may 

run  into  problems  if  either  p0  or  q0  turns  out  to  be  zero.  For  example,  when  J3  =  1  and  po  =  -1,  qo 
=  1,  then  q0  =  (0,  2(1  +/3po),  0)  =  (0,  0,  0).  If  we  define  the  unit  normal  in  terms  of  these  two  vec¬ 
tors,  then  in  this  example,  we  will  run  into  the  case  of  dividing  by  zero.  In  such  a  case,  and  since 
the  unit  normally  does  not  change  on  a  flat  BQ  we  can  do  either  of  two  things:  select  two  non- 
collinear  vectors  on  the  BQ  to  define  the  unit  normal  or,  even  better,  ignore  it  altogether  since  it 
only  appears  in  B  in  (7.4)  and,  there,  it  is  dotted  to  the  vector  rpq  which  lies  on  the  BQ's  plane; 
thus,  this  product  is  zero  and  can  be  omitted. 

In  Tables  9.3  and  9.4,  we  present  results  for  the  BQ  of  (9.20)  and  for  three  different  values 
of  p.  In  Tables  9.3,  we  take  the  singular  point  to  be  the  origin.  We  see  that  the  results  from  the 
two  quadratures  agree  to  at  least  14SD.  We  also  note  that  the  results  of  Tables  9.3.1  and  9.3.3  are 
identical.  This  is  because  the  BQ  (a  triangle  in  this  case)  is  the  reflection  of  the  BQ  in  the  second 
case  about  the  y-axis,  as  shown  in  Figs.  9.3.1  and  9.3.2.  If  we  compare  the  values  in  Table  9.1.3 
to  those  in  Table  9.3.1,  we  see  a  smooth  transition  between  the  values  for  ft  =  0.999  and  those  for 
p  =  1.  We  show  the  BQ  for  Table  9.3.2  in  Figure  9.3.3.  It  comprises  two  triangles  with  a  com¬ 
mon  vertex. 

Table  9.3.1:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  p0  =  q0  =  0 
and  P=  1. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.74737108607702 

0.031 

15 

2.74737108607702 

0.031 

15 

1.00E-08 

2.85044223392429 

0.046 

14 

2.85044223392428 

0.047 

14 

1.00E-20 

2.85044224439626 

0.000 

15 

2.85044224439626 

0.016 

15 

0.00E+00 

2.85044224439626 

0.031 

15 

2.85044224439626 

0.000 

15 

Table  9.3.2:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  qo  =  0 
and  P=  2. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.73594640126082 

0.016 

15 

2.73594640126082 

0.016 

15 

1.00E-08 

2.83994915267823 

0.016 

15 

2.83994915267823 

0.031 

15 

1.00E-20 

2.83994916315021 

0.000 

15 

2.83994916315021 

0.015 

15 

0.00E+00 

2.83994916315021 

0.016 

15 

2.83994916315021 

0.031 

15 
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Table  9.3.3:  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  qo  =  0 
and  /?  =  - 1 . 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

2.74737108607702 

0.016 

15 

2.74737108607702 

0.015 

15 

1.00E-08 

2.85044223392429 

0.031 

14 

2.85044223392428 

0.031 

14 

1.00E-20 

2.85044224439626 

0.031 

15 

2.85044224439626 

0.000 

15 

0.00E+00 

2.85044224439626 

0.000 

15 

2.85044224439626 

0.016 

15 

y 


Figure  9.3.1:  Graph  of  the  BQ  of  (9.19)  with  /?=  1  (p'  =  -1.0). 
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y 


Figure  9.3.2:  Graph  of  the  BQ  of  (9.19)  with  /?  =  -1  (p'=  1.0). 


-3-2-10  1  2  3 


Figure  9.3.3:  Graph  of  the  BQ  of  (9.19)  with  /?=  2  (//=  -0.5). 
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In  Tables  9.4,  we  take  the  singular  point  to  be  the  point  ( po  =  -1,  qo  =  1).  In  Table  9.4.1,  we 
present  the  results  for  J3=  1.  In  this  case,  the  singular  point  is  the  point  (-3,  0,  0)  in  Figure  9.3.1. 
Also,  we  have  that  qo  =  (0,  0,  0),  as  we  mentioned  above.  In  this  case,  we  have  from  (7.4)  that 

A  =  h2  ,  B  =  0,  C  =  ^j 9  +  4(l  +  v)2  .  (9.23) 

When  we  substitute  these  values  into  the  logarithm's  argument  in  (7.9),  we  get 


7°  l 

1=  f  — Tln 


[2C(v)]~  +h2  +  2  C(v) 


1*1 


dv 


(9.24) 


an  expression  that  is  not  computable  when  h  =  0.  If  we  compare  this  table  to  Table  9.2.3,  we  see 
that  the  transition  for  h  =  0.1  is  smooth.  As  h  gets  smaller,  however,  we  see  a  radical  change  in 
the  integration  outcome.  We  note  that  in  Table  9.2.3,  the  value  of  qo  is  (0,  0.002,  0),  a  small  val¬ 
ue  but  sufficient  to  stabilize  the  integral. 

Table  9.4.1.  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  -1,  qo  =  1 
and  /?=  1. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

3.03215129971745 

0.015 

14 

3.03215129971744 

0.016 

14 

1.00E-08 

13.1082615780428 

0.000 

15 

13.1082615780428 

0.015 

15 

1.00E-20 

30.3816595133087 

0.031 

15 

30.3816595133087 

0.016 

15 

0.00E+00 

Failed  to  compute 

N/A 

N/A 

Failed  to  compute 

N/A 

N/A 

Table  9.4.2.  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  -1,  qo  =  1 
and  /?  =  2. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

1.27899151358520 

0.000 

14 

1.27899151358519 

0.000 

14 

1.00E-08 

1.29041253658337 

0.016 

14 

1.29041253658336 

0.031 

14 

1.00E-20 

1.29041253765587 

0.000 

15 

1.29041253765587 

0.015 

15 

0.00E+00 

1.29041253765587 

0.015 

15 

1.29041253765587 

0.031 

15 
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Table  9.4.3.  The  value  of  the  integral  in  (7.9)  for  the  BQ  of  (9.19)  and  for  po  =  -1,  r/o  =  1 
and  /?=-!. 


h 

GKQ 

DEQ 

Value 

Time 

SD 

Value 

Time 

SD 

1.00E-01 

0.936274735416838 

0.015 

14 

0.936274735416837 

0.016 

14 

1.00E-08 

0.944587582282624 

0.015 

14 

0.944587582282623 

0.016 

14 

1.00E-20 

0.944587583101619 

0.016 

14 

0.944587583101618 

0.000 

14 

0.00E+00 

0.944587583101619 

0.000 

14 

0.944587583101618 

0.016 

14 

In  concluding  this  section,  we  stress  that  we  should  avoid  irregular-  BQs  in  practice.  For  flat 
but  regular  BQs,  we  can  use  (7.9)  without  incurring  a  loss  of  SD. 
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SECTION  10:  SUMMARY  AND  CONCLUSIONS 

We  have  presented  a  new  approach  for  computing  the  simple  layer  potential  of  polynomial 
density  over  a  BQ.  The  analysis  focuses  on  the  case  when  the  OP  is  near  or  on  the  BQ.  Our  ap¬ 
proach  differs  from  existing  ones  in  that  it  does  not  use  a  transformation  of  coordinates  nor  does 
it  employ  an  auxiliary  plane;  instead  it  uses  rectangular  coordinates  and  integrates  analytically 
along  one  of  them  and  numerically  along  the  other.  We  have  demonstrated  that,  with  carefully 
designed  quadratures,  this  approach  can  yield  DP  for  any  regular  BQ  and  any  position  of  the  OP. 

In  Section  2,  we  define  the  BQ  mathematically  and  provide  its  geometric  properties.  We  al¬ 
so  make  the  distinction  between  a  regular  and  an  irregular  BQ;  moreover,  we  show  that  an  irreg¬ 
ular  BQ  is  necessarily  flat  and  we  provide  examples.  In  Section  3,  we  describe  three  BQs  that  we 
use  for  testing  our  approach.  They  range  from  a  flat  one  to  a  very  elongated  three-dimensional 
one. 


In  Section  4,  we  discuss  in  detail  the  reasons  for  the  need  of  a  new  approach  to  the  calcula¬ 
tion  of  the  simple  layer  potential.  We  demonstrate  through  analysis  and  examples  that,  under  the 
existing  approach  of  the  tangent-plane  approximation,  the  limit  of  the  integrand  does  not  exist  at 
the  singularity  point,  and  that  the  integrand  may  undergo  rapid  changes  in  the  neighborhood  of 
this  point.  We  also  demonstrate  that  the  use  of  polar  coordinates  works  well  when  the  OP  is  on 
the  BQ  but  not  when  it  is  near  the  BQ. 

In  Section  5,  we  introduce  our  approach.  We  split  the  integrand  into  two  parts:  one  that  has 
no  singularities  of  any  kind,  and  one  that  does  not  oscillate  but  contains  the  singularity.  The  inte¬ 
gral  of  the  formal  can  be  evaluated  numerically  to  good  precision.  For  the  latter  we  show  that, 
through  repeated  integration  by  parts,  we  can  reduce  its  integral  to  a  number  of  single  integrals 
with  at  least  continuous  integrands  and  one  that  is  basically  a  simple  layer  (or  Newtonian  poten¬ 
tial)  of  density  one. 

In  Section  6,  we  consider  the  simple  layer  of  density  one  and  integrate  it  with  respect  to  one 
of  the  two  parameters  that  define  the  BQ.  Through  analysis  and  examples,  we  show  that  the  re¬ 
maining  integral  (with  respect  to  the  remaining  parameter  and  which  can  be  evaluated  only  nu¬ 
merically)  loses  precision  as  the  OP  approaches  the  BQ.  In  Section7,  we  study  the  cause  of  this 
behavior  and  find  that  it  is  due  to  two  terms  in  the  integrand  that  subtract  one  another.  When  the 
integration  variable  is  equal  to  zero,  and  as  the  OP  approaches  the  BQ,  the  two  terms  become 
almost  equal  and  their  difference  appears  as  zero  while  it  is  not,  unless  the  OP  is  on  the  BQ.  We 
correct  for  this  by  rewriting  the  expression  for  the  integrand  in  such  a  way  that  the  difference  of 
the  offending  terms  is  multiplied  by  a  term  that  becomes  zero  when  the  integration  variable  is  ze¬ 
ro.  Using  examples,  we  show  that  the  new  expression  is  stable  and  that  we  can  obtain  15  SD  of 
precision  even  when  the  OP  lies  on  the  BQ.  In  Section  8,  we  continue  testing  our  method  using 
the  third  BQ  of  Section  3.  For  all  three  test  points  we  get  a  minimum  of  14  SD.  We  also  use  inte¬ 
gration  by  parts  to  produce  a  new  formula  that  is  guaranteed  to  be  stable  when  the  OP  is  on  the 
BQ.  As  we  point  out  there,  all  numerical  experiments  with  the  general  formula  yield  a  minimum 
of  14  SD;  moreover,  even  though  the  OP  lies  on  the  BQ,  we  can  move  it  a  small  distance  off  it 
without  loss  of  precision.  Thus,  we  may  not  need  to  use  this  special  formula,  a  formula  that  is 
more  computationally  intensive  than  the  general  one. 
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In  Section  9,  we  prove  a  claim  we  made  in  earlier  sections  and  show  that  a  certain  term  can 
be  zero  only  if  the  BQ  is  irregular.  We  also  give  additional  examples  of  irregular  BQs  and  per¬ 
form  a  sensitivity  study  to  determine  how  precision  is  affected  as  we  pass  from  a  regular  to  an  ir¬ 
regular  BQ.  As  we  point  out  there,  irregular  BQs  should  be  avoided  in  practice. 

In  terms  of  additional  investigations,  we  may  wish  to  perform  a  more  detailed  validation  of 
our  formulas.  We  believe,  however,  that  they  will  withstand  any  test,  no  matter  how  severe.  We 
may  also  perform  calculations  for  several  OPs  lying  on  the  same  BQ  to  determine  the  resulting 
surface  and  its  smoothness.  This  is  necessary  for  determining  whether  the  second  integration,  the 
one  with  respect  to  the  OP  coordinates  can  be  performed  numerically  to  high  precision.  Alterna¬ 
tively,  we  may  wish  to  use  these  points  to  interpolate  and  use  the  resulting  expressions  for  an  an¬ 
alytic  integration  with  respect  to  the  OP  coordinates. 
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APPENDIX  A:  Irregular  Bilinear  Quadrilateral 

In  this  appendix,  we  prove  the  oral  arguments  that  follow  (2.8).  We  resolve  the  vector  rpq 
with  a  component  in  the  plane  formed  by  the  vectors  rp  and  rq  and  a  component  perpendicular  to 
that  plane 

Xpq  =  arp  +  PXq  +  7  (rp  X  Ty )  •  (A- 1) 

In  place  of  (2.8),  we  can  then  write 


dr(p,q)  Jr (p,g)=r  xr  + 
dp  dq  p  q 


avp  +  p rq  +  y  (rp  x  rq )]  x  (r qq  -  r pp) 

=  (\  +  aq  +  pp)rpxrq+r (rp  xrjx (r qq - r pP) 

=  (l  +  aq  +  pp)(rpxrq)  +  r\rp  p(rp -rq)~ q\rq\~  +rq  q(rp-rq)-p\r 


(A.2) 


This  is  a  vector  for  a  component  in  the  plane  of  the  vectors  rp  and  rq,  and  one  perpendicular  to 
this  plane.  For  this  vector  to  be  equal  to  zero,  all  three  components  must  be  zero 


1  +  aq  +  ftp 

=  0 

r  / 

2" 

=  0 

1 

y 


(A.3) 


If  we  assume  that  y  is  different  from  zero,  then  the  last  two  equations  yield  a  determinant  differ¬ 
ent  from  zero  since  we  have  assumed  that  rp  and  rq  are  not  co-linear.  Thus,  the  system  of  two 
equations  has  only  the  trivial  solution,  i.e.,  p  and  q  are  both  equal  to  zero.  For  these  values,  how¬ 
ever,  the  first  equation  in  (A.3)  is  not  satisfied.  This  leads  us  to  the  conclusion  that,  for  the  two 
equations  to  be  equal  to  zero,  y  must  be  equal  to  zero.  From  (A.l),  we  conclude  that,  for  a  BQ  to 
be  irregular,  the  three  vectors  describing  it  must  lie  on  the  same  plane.  If  we  assume  that  all  three 
are  non-zero  vectors,  we  conclude  that  a  necessary  condition  for  an  irregular  triangle  is  that 

Vpq  Xr9)=0’  (A‘4) 

In  case  rp  and  rq  are  colinear,  the  first  tenn  on  the  right  of  (2.8)  is  zero  and  the  second  be¬ 
comes  zero  at  the  origin;  thus,  we  have  an  irregular  BQ  that  lies  on  the  plane  defined  by  rp  (or  rq) 
and  rpq. 
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APPENDIX  B:  Continuation  of  Analysis  of  Section  4 

In  Section  4,  we  analyzed  the  inner  integral  of  the  first  BQ  of  Section  3.  We  found  that  the 
first  integrand  in  (4.2)  is  not  as  smooth  as  stated  in  [4].  Here,  we  perform  the  same  analysis  for 
the  remaining  two  BQs  of  Section  3. 

The  second  BQ  is  shown  in  Figure  3.2.  The  position  vector  is  given  by  (3.6)  so  that 

r*  =  (3, 0, 0)  \  +  (0, 2, 0)  \  +  (0, 0,  l)  (B .  1 ) 

4  4  16 


and,  hence, 

r  («,  v)  -  r*  =  (3, 0, 0)  u  +  (0, 2, 0)  v  +  (0, 0,  l) 
The  linear  part  of  this  is 


r0(n,v)  -r*  =  (3,0,0)u  +  (0, 2,0)v  +  (0,0,1)^-  (u  +  v)  . 
From  the  last  two  we  get  that 


7  O 

(  O 

1  ' 

u- 1 — 

V  H - 

— 

LI  4  J 

V  4  J 

16 

R2(u,v)=  (3u)2  +(2  v)2  + 


'(  n 

c  n 

1  " 

+ 

V  H - 

— 

l  4  J 

v  4 ) 

!6 

and 


Ro,2(u’v)  =  J(3«)2+(2v)2  + 


^  u  + 


V  *  J 


(B.2) 


(B.3) 


(B.4) 


(B.5) 


so  that 


h2  (m,v)  =  ■ 


(  n 

/ 

r  n 

v+- 

l  4  J 

1 

l  4  J 

1 

v  4/ 


(B.6) 


|(3m)2  +  (2v)2  + 


Y  n 

(  O 

1  1 " 

+ 

u  +  — 

V  H - 

\  4 ) 

v  4  J 

16 

(3 u)  +(2v)  + 


u  +  v 


We  have  plotted  this  expression  in  Mathematica  [11].  Figure  B.l  shows  the  graph  of  (B.6),  a 
graph  that  is  anything  but  smooth  near  the  origin.  In  Figs.  B.2  and  B.3,  we  zoom  toward  the 
origin  and  take  cuts  along  both  axes.  We  see  that,  on  at  least  one  side  of  an  axis,  we  have  a  fast¬ 
changing  graph  that  is  difficult  to  match  numerically. 
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Figure  B.l:  Graph  of  the  surface  (B.6)  for  1  -  0 ,m  =  0 ;  Second  BQ. 


Figure  B.2.  Cuts  of  Fig.  B.l  along  v  =  -0.05 (blue),  v  =  0  (invisible  but  coinciding  with  the  u- 
axis)  and  v  =  0.05  (olive).  Second  BQ. 
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Figure  B.3.  Cuts  of  Fig.  B.l  along  u  = -0.05(blue),  u  =  0  (invisible  but  coinciding  with  the  v- 
axis)  and  u  =  0.05  (olive).  Second  BQ. 

In  Figs.  B.4  and  B.5,  we  display  the  behavior  of  (B.6)  for  /  =  l,m  =  0  while,  in  Figs.  B.6  and 
B.7,  we  display  the  graph  of  the  same  function  for  /  =0,m  =  l.  In  Figs.  B.8  -  B.ll,  we  display 
the  behavior  for  /  =  l,m  =  1.  What  we  observe  is  very  similar  to  what  we  discussed  in  Section  4 
for  the  first  BQ  of  Section  3.  We  will  not  repeat  those  comments  here. 
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Figure  B.4:  Graph  of  (B.6)  for  l  =  l,m  =  0 ;  Second  BQ. 
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Figure  B.5:  Graph  of  the  function  h2  (w,0) ,  as  defined  in  (B.6),  for  /  =  1 ,  m  =  0 ;  Second  BQ. 
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Figure  B.6:  Graph  of  (B.6)  for  /  =  0,m  =  1 ;  Second  BQ. 
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Figure  B.7:  Graph  of  the  function  h2  (0,  v)  ,  as  defined  in  (B.6),  for  /  =  0,m  =  1 ;  Second  BQ. 


Figure  B.8:  Graph  of  (B.6)  for  /  =  l,m  =  1;  Second  BQ. 
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Figure  B.9:  Graph  of  the  function  h2  (w,0) ,  as  defined  in  (B.6),  for  /  =  l,m  =  1 ;  Second  BQ. 
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Figure  B.10:  Graph  of  the  function  h2  (0,v)  ,  as  defined  in  (B.6),  for  /  =  l,m  -  1 ;  Second  BQ. 


Figure  B.ll:  Graph  of  the  function  h2  (u,u),  as  defined  in  (B.6),  for  /  =  l,m  =  1;  Second  BQ. 
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We  proceed  next  to  the  last  BQ  of  Section  3.  Its  position  vector  is  given  by  (3.9).  From  it 


r  *  =  (3, 0, 0)  ^  +  (0, 2, 0)  ^  +  (0, 0 ,10)  ^ 


and,  hence, 

r (p,  q)  - r*  =  (3, 0, 0) u  +  (0, 2, 0)  v  +  (0, 0, 10) 
The  linear  part  of  this  is 


\(  0 

(  n 

1  ' 

V  H - 

— 

LI  4  J 

v  4  J 

16 

1 


r0  (u,  v)  -  r  *  =  (3, 0, 0)  w  +  (0, 2, 0)  v  +  (0, 0, 1 0)  —  (w  +  v) 
From  the  last  two  we  get  that 


R3  (w,v)  =  |r(p,g)  -r  *|  =  (3u)2  +  (2 vf  +100 


\(  O 
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|  1  " 

v  + 

— 

LI  4  J 
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4  J 

~16_ 

“|2 


(B.7) 


(B.8) 


(B.9) 


(B.10) 


and 


/?30(M,v)  =  |r0(p,^)-r*|  =  J(3m)2+(2v)2+  -  (u  +  vf 

\2J 


(B.  11) 


so  that 


h3  (w,v)  = 


f  l') 
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f  1) 

m 

rn 

— 

l  4  J 

1 

l  4  J 

U  ) 

l+m 


(3m)2+(2v)2+100 


(  n 

l  4  J 


(  0 

1  1 

T  f, 

V  H - 

l  4  J 

!6  J 

V 

■  (B.12) 


(3u)2+(2v)2  + 


—  (u  +  v) 
2V  ’ 


We  use  Mathematica  [11]  to  produce  graphs.  In  Figure  B.12  we  show  the  graph  of  (B.6),  a 
graph  that  is  anything  but  smooth  near  the  origin.  In  Figs.  B.13  and  B.14,  we  zoom  toward  the 
origin  and  take  cuts  along  both  axes.  We  see  that,  on  either  side  of  an  axis,  we  have  a  doublet¬ 
like  behavior  that  is  difficult  to  match  numerically. 
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Figure  B. 12:  Graph  of  the  surface  (B.  12)  for  / -0,m  =  0 ;  Third  BQ. 


Figure  B.13.  Cuts  of  Fig.  B.12  along  v  =  -0.05  (blue),  v  =  0  (invisible  but  coinciding  with  the  u- 
axis)  and  v  =  0.05  (olive).  Third  BQ. 


88 


APPENDIX  B 


NAWCADPAX/TR-20 15/214 


Figure  B. 14.  Cuts  of  Fig.  B.  12  along  u  =  -0.05 (blue),  u  =  0  (invisible  but  coinciding  with  the  v- 
axis)  and  u  =  0.05  olive.  Third  BQ. 


In  Figs.  B.15  and  B.16,  we  display  the  behavior  of  (B.12)  for  I  =  1 , m  =  0  while,  in  Figs. 
B.17  and  B.18,  we  display  the  graph  of  the  same  function  for  /  =  0,/m  =  1 .  In  Figs.  B.19  -  B.22, 
we  display  the  behavior  for  /  =  l,m  =  1 .  What  we  observe  is  very  similar  to  what  we  discussed  in 
Section  4  for  the  first  BQ  of  Section  3,  except  for  the  last  graph,  where  we  display  the  behavior 
along  the  diagonal.  We  verify  below  that  the  graph  displays  the  correct  behavior.  From  (B.12) 
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which  explains  the  non-linearity  of  the  graph  in  Figure  B.22.  Taking  the  limit,  we  have 


This  number  is  in  agreement  with  the  one  we  get  from  Figure  B.22. 

The  bottom  line  here  is  that  the  integrand  of  the  first  integral  in  (4.2)  is  not  sufficiently 
smooth  and  that  a  numerical  evaluation  of  it  would  require  very  many  grid  points. 
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Figure  B.15:  Graph  of  (B.12)  for  l  =  \,m  =  0;  Third  BQ. 


Figure  B.16:  Graph  of  the  function  /z3  [u, 0)  ,  as  defined  in  (B.12),  for  /  =  l,m  -  0 ;  Third  BQ. 
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V 


Figure  B.17:  Graph  of  (B.12)  for  /  =  0 ,m  =  1 ;  Third  BQ. 
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Figure  B.18:  Graph  of  the  function  h3  (0,  v) ,  as  defined  in  (B.12),  for  /  =  0,m  =  1 ;  Third  BQ. 
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Figure  B.19:  Graph  of  (B.12)  for  /  =  1  ,m  =  1;  Third  BQ. 
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Figure  B.20:  Graph  of  the  function  /z3  [u, 0)  ,  as  defined  in  (B.12),  for  /  =  l,m  =  1 ;  Third  BQ. 
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Figure  B.21:  Graph  of  the  function  /?.,  (0,  v) ,  as  defined  in  (B.12),  for  /  =  l,m  =  1 ;  Third  BQ. 


Figure  B.22:  Graph  of  the  function  h 3  (m,«)  ,  as  defined  in  (B.12),  for  /  =  l,m  =  1 ;  Third  BQ. 


94 


APPENDIX  B 


NAWCADPAX/TR-20 15/214 


DISTRIBUTION: 


NAVAIRSYSCOM  (AIR-4.5. 5/Douglas  McLaughlin),  Bldg.  2187,  Room  3201  (1) 

48110  Shaw  Road,  Patuxent  River,  MD  20670-1906 
WIPL-D  d.o.o.  (Prof.  Branko  Kolundzija),  Gandijeva  7  apt.  32,  (1) 

11073  Belgrade,  Serbia  (branko.kolundzija@wipl-d.com) 

University  of  Houston  (Prof.  Donald  R.  Wilton),  E.E.  Dept.  (1) 

N  308  Engineering  Building  1,  Houston,  TX  77204-4005 
NAVAIRSYSCOM  (AIR-5. IV),  Bldg.  304,  Room  106A  (1) 

22541  Millstone  Road,  Patuxent  River,  MD  20670-1606 
NAVTESTWINGLANT  (55TW01A),  Bldg.  304,  Room  200  (1) 

22541  Millstone  Road,  Patuxent  River,  MD  20670-1606 
NAVAIRSYSCOM  (AIR-5.1),  Bldg.  304,  Room  100  (1) 

22541  Millstone  Road,  Patuxent  River,  MD  20670-1606 
NAVAIRSYSCOM  (Air-4.0T),  Bldg.  407,  Room  116  (1) 

22269  Cedar  Point  Road,  Patuxent  River,  MD  20670-1120 
DTIC  (1) 

Suite  0944,  8725  John  J.  Kingman  Road,  Ft.  Belvoir,  VA  22060-6218 


95 


